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STATISTICAL  PREDICTION  OF  LAMINAR-TURBULENT  TRANSITION 


ROBERT  RUBINSTEIN"  AND  MEELAN  CHOUDHARF 

Abstract.  Stochastic  versions  of  stability  equations  are  considered  as  a  means  to  develop  integrated 
models  of  transition  and  turbulence.  Two  types  of  stochastic  models  are  considered:  probability  density 
function  evolution  equations  for  stability  mode  amplitudes,  and  Langevin  models  based  on  representative 
stability  theories  including  the  resonant  triad  model  and  the  parabolized  stability  equations.  The  first  type 
of  model  can  describe  the  effect  of  initial  phase  differences  among  disturbance  modes  on  transition  location. 
The  second  type  of  model  describes  the  growth  of  random  disturbances  as  transition  proceeds  and  provides 
a  natural  framework  in  which  to  couple  transition  and  turbulence  models.  Coupling  of  parabolized  stability 
equations  with  either  subgrid  stress  models  or  with  conventional  turbulence  models  is  also  discussed  as  an 
alternative  route  to  achieve  the  goal  of  integrated  turbulence  and  transition  modeling. 

Key  words,  turbulence,  transition,  parabolized  stability  equations,  resonant  triads 

Subject  classification.  Fluid  Mechanics 

1.  Introduction.  Recent  developments  in  computational  fluid  dynamics  flows  have  revealed  the  need 
for  integrated  modeling  of  turbulence  and  transition.  In  many  external  aerodynamic  flows,  such  as  the  flow 
past  high-lift  airfoil  configurations,  low-pressure  turbine  blades,  and  the  Mars  flyer  design  proposed  by  NASA 
(C.  L.  Streett,  private  communication),  accurate  prediction  of  crucial  bulk  parameters  like  the  lift  and  drag 
coefficients  requires  that  flow  computation  be  transition  sensitized.  In  high-lift  applications  in  particular,  the 
laminar  and  turbulent  regions  are  tightly  coupled  through  the  dependence  of  separation  and  reattachment 
characteristics  on  transition  location;  the  predicted  flow  can  vary  from  fully  attached  to  massively  separated 
depending  on  the  transition  location  imposed  on  the  how  solver  [25].  Recent  computations  [38]  have  shown 
that  whereas  most  turbulence  models  can  predict  many  external  hows  with  adequate  accuracy  provided  a 
transition  location  is  specified  in  advance,  no  single  turbulence  model  can  predict  this  location  in  all  cases. 

A  basic  obstacle  to  the  integration  of  turbulence  and  transition  modeling  lies  in  the  quite  distinct 
viewpoints  of  previous  research  in  these  helds.  Transition  studies  typically  formulate  deterministic  equations 
like  the  Orr-Sommerfeld  equation  in  a  wide  variety  of  different  problems,  whereas  turbulence  studies  apply 
statistical  methods  to  describe  the  features  common  to  all  turbulent  hows.  This  distinction  is  best  illustrated 
by  comparing  the  simplest  turbulence  model,  the  mixing-length  model  [49]  and  the  simplest  transition 
prediction  scheme,  the  eN  method  [5],  [41].  The  eA  method  is  a  deterministic  procedure  based  on  the  linear 
amplification  of  the  most  unstable  mode,  whereas  the  mixing-length  model  appeals  to  a  statistical  mechanics 
analogy  to  describe  the  effects  of  turbulence  on  the  mean  how. 

Although  these  divergent  viewpoints  arise  from  the  need  to  describe  radically  different  how  physics, 
recent  developments  have  tended  to  bring  these  helds  closer  together;  thus,  recent  investigations  of  bypass 
transition  [33]  analyze  random  initial  disturbances,  implicitly  posing  the  transition  problem  statistically.  In 
any  case,  a  synthesis  of  turbulence  and  transition  modeling  is  demanded  by  the  need  for  accurate,  seamless 
computation  of  transitional  hows  throughout  the  entire  transition  process,  from  the  initial  disturbances  to 
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fully  developed  turbulence. 

One  approach  to  integrated  modeling  of  turbulence  and  transition  is  to  apply  low  Reynolds  number 
turbulence  models  directly  to  transitional  flows  [1],  [7],  [11];  this  approach  treats  transition  as  a  special 
kind  of  ‘non-equilibrium’  turbulence.  Another  general  approach  is  to  sensitize  existing  turbulence  models 
to  transition,  often  by  introducing  a  new  field  quantity  like  intermittency  [45],  [26],  [40].  An  entirely  new 
approach  which  treats  transition  beginning  from  turbulence  modeling  is  given  in  [46].  Certainly,  models  of 
this  type  will  play  a  valuable  role  in  the  aerodynamic  design  process.  However,  due  to  the  limited  transition 
physics  embedded  in  such  models,  typically  through  curve  fits  and/or  ad  hoc  calibration  against  a  relatively 
sparse  measurement  database,  they  are  not  likely  to  explain  potential  anamolies  in  measured  aerodynamic 
performance  or  yield  reliable  answers  when  extrapolation  to  completely  new  flow  situations  (e.g.,  the  Mars 
flyer  application)  is  required. 

In  broad  terms,  the  present  work  is  motivated  by  the  desire  to  retain  as  much  of  the  transition  physics 
revealed  by  recent  advances  in  understanding  of  the  transition  process  reviewed,  for  example  by  [16],  [22], 
[14],  and  [27],  yet  allow  for  integrated  aerodynamic  predictions  in  a  realistic  setting.  To  that  end,  the  current 
investigation  begins  with  transition  and  has  the  ultimate  goal  of  incorporating  turbulence  prediction  in  a 
seamless  fashion.  Given  this  overall  goal,  two  main  areas  present  themselves  as  candidates  for  further  work. 
The  first  of  these  areas  concerns  the  prediction  of  transition  in  isolation  from  the  turbulent  flow  computation, 
whereas  the  second  area  of  need  is  related  to  the  mechanics  of  coupling  the  transition  prediction  with  the 
computation  of  the  turbulent  flowfield. 

First,  let  us  consider  the  prediction  of  transition  itself.  As  noted  by  Reshotko  [34],  laminar-turbulent 
transition  in  a  convectively  unstable  boundary  layer  corresponds  to  the  forced  response  of  a  rather  complex 
oscillator,  namely  the  laminar  basic  state.  Under  the  relatively  benign  disturbance  environments  that  are 
typical  of  external  aerodynamic  flight,  the  above  response  may  be  decomposed  into  three  broad  stages: 
excitation  of  linearly  unstable  eigenrriodes  of  the  oscillator  due  to  forcing  from  the  external  disturbances  (i.e., 
the  process  of  receptivity),  initially  exponential  growth  of  these  eigenmodes  as  they  propagate  downstream, 
and  the  variety  of  nonlinear  interactions  which  ensue  after  the  eigenmodes  have  achieved  sufficiently  large 
amplitudes  and  which  eventually  lead  to  turbulence.  Any  attempt  to  model  the  transition  process  on  a 
rational  basis  must  account  for  all  three  stages  of  the  ‘oscillator’  response. 

The  current  engineering  methods  for  transition  prediction  typically  involve  the  linear  propagation  phase 
alone,  with  little  or  no  account  of  receptivity  (which  initiates  transition)  or  the  nonlinear  interactions  (which 
actually  cause  the  onset  of  transition).  In  spite  of  their  success  in  predicting  transition  in  a  wide  range  of 
flows,  these  JY-factor  methods  lack  the  physical  basis  necessary  to  explain  the  effects  of  surface  roughness 
and  free-stream  unsteadiness,  which  may  be  particularly  significant  in  flows  such  as  swept  wing  boundary 
layers  and  high-lift  flows.  Thus,  as  an  essential  aspect  of  integrated  modeling  of  transition  and  turbulence, 
it  is  necessary  to  develop  rational  prediction  methods  that  couple  the  various  stages  of  transition  itself. 
A  preliminary  attempt  along  these  lines,  which  focused  on  a  controlled  disturbance  environment  involving 
sinusoidal  surface  waviness  and  a  monochromatic  acoustic  disturbance,  was  presented  in  [3].  However, 
considerable  advancement  is  necessary  to  enable  similar  predictions  in  a  realistic  setting. 

In  view  of  the  complexity  of  the  forcing  fields  encountered  in  practical  applications,  integrated  transition 
prediction  needs  to  be  performed  within  a  stochastic  framework  that  is  capable  of  analyzing  the  effects 
of  randomness  and  uncertainties  in  the  flow  environment.  The  concept  of  stochastic  transition  prediction 
represents  a  major  theme  behind  the  work  reported  in  this  paper.  The  objective  here  is  not  so  much  to 
develop  and  present  a  well-tested  recipe  towards  this  goal,  but  to  outline  the  various  issues  introduced  by 
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stochasticity,  explore  potential  ways  in  which  they  might  be  addressed,  and  examine  the  potential  merits 
and  demerits  of  each  approach.  Although  mainly  aimed  at  integrated  transition  and  turbulence  modeling, 
we  note  that  such  work  would  also  be  relevant  to  developing  stochastic  control  schemes  for  laminar-turbulent 
transition. 

In  principle,  statistical  consideration  of  the  forced-oscillator  response  should  include  the  uncertainties 
associated  with  both  oscillator  characteristics  (i.e.,  undisturbed  laminar  state)  and  the  forcing  environment. 
Of  these  two,  the  stochastic  character  of  the  forcing  is  perhaps  more  important,  and  is  what  is  emphasized 
here.  In  general  terms,  we  consider  how  deterministic  models  for  a  given  stage  of  the  transition  process 
(which  may  be  considered  as  well-established  at  this  point)  may  be  extended  to  a  stochastic  setting  by 
treating  the  input  conditions  as  random  rather  than  fully  specified. 

In  this  paper,  we  assume  that  the  dominant  receptivity  occurs  through  linear  mechanisms  and  hence  the 
overall  problem  of  stochastic  disturbance  evolution  in  a  boundary  layer  may  be  decomposed  into  two  separate 
problems:  the  receptivity  problem  that  predicts  the  statistics  of  initial  disturbances  (i.e.,  instability  modes) 
as  functions  of  input  forcing,  and  the  propagation  stage  which  predicts  how  these  statistical  characteristics 
evolve  as  the  instability  waves  amplify  and  eventually  interact  through  nonlinear  mechanisms.  Given  the 
system  parameters  (i.e.,  a  plant  model),  the  response  of  a  linear,  deterministic  system  to  stochastic  input 
can  be  predicted  rather  simply  since  the  second  moments  of  the  response  are  directly  proportional  to  those 
of  the  input,  with  the  scaling  factor  begin  the  known  (deterministic)  transfer  function.  Thus,  it  is  relatively 
straightforward  to  extend  the  available  models  for  receptivity  and  linear  growth  to  a  stochastic  setting. 
Accordingly,  we  focus  our  attention  on  stochastic  aspects  of  the  nonlinear  propagation  phase,  for  which 
there  is  no  unique  transfer  function  as  the  latter  now  becomes  a  function  of  the  input  parameters  (i.e.,  initial 
amplitude  spectrum). 

As  transition  progresses,  a  sequence  of  secondary  and  higher  order  instabilities  are  excited  leading  to  the 
emergence  of  smaller  scales  of  motion;  with  continued  excitation  of  small  scales,  memory  of  the  details  of  the 
original  state  fades,  and  gradually  a  statistically  steady  state  of  turbulence  is  reached.  Two  points  are  worth 
noting  in  this  context.  First,  as  the  system  becomes  chaotic  during  the  laminar  breakdown  stage,  stochastic 
modeling  is  necessary  even  if  there  were  no  uncertainties  in  the  forcing  environment.  Second,  stochastic 
formulation  of  the  transition  models  provides  a  natural  link  between  transition  theory  and  turbulence  mod¬ 
eling.  Of  course,  passing  from  a  deterministic  to  a  stochastic  transition  analysis  is  not  trivial,  particularly 
in  view  of  nonlinearity  of  the  problem.  The  details  of  the  path  to  turbulence,  in  which  different  modes  are 
significant  during  different  phases  in  the  evolution,  is  highly  case  dependent.  Therefore,  a  direct  simulation 
of  the  stochastic  problem,  which  would  require  simulating  a  large  number  of  these  paths,  is  impractical.  This 
problem  is  much  more  severe  than  the  problem  of  simulating  fully  developed  turbulence,  in  which  ergodicity 
helps  permit  the  evaluation  of  statistics.  The  present  paper  represents  but  a  first  step  to  a  statistical  model 
of  transition. 

An  outline  of  the  paper  is  as  follows.  To  begin,  we  briefly  review  the  theory  of  receptivity,  since  in 
many  external  flows,  receptivity  analysis  will  determine  the  initial  disturbance  spectrum.  Following  [8],  we 
emphasize  the  stochastic  aspects  of  the  problem,  namely  the  connection  between  the  random  free-stream 
disturbances  and  roughness  distributions  and  the  probability  density  function  of  the  phase  and  amplitude  of 
the  generated  Tollmien-Schlichting  waves. 

The  subsequent  evolution  of  the  generated  instability  waves  (particularly,  during  the  nonlinear  stage)  is 
examined  next.  On  a  determinstic  level,  the  disturbance  evolution  problem  can  be  addressed  using  models 
of  varying  degrees  of  sophistication  and  accuracy  including  Craik’s  model  for  resonant  triads  of  Tollmien- 
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Schlichting  waves  [12],  rational  high  Reynolds  number  asymptotics  [14],  parabolized  stability  equations  [4], 
large  eddy  simulation  [18],  and  direct  numerical  simulation  [33], [51].  To  achieve  the  goal  of  integrated 
transition  and  turbulence  modeling,  we  start  by  applying  these  models  in  a  stochastic  setting,  treating  the 
initial  disturbance  spectrum  as  random  rather  than  deterministic.  We  begin  by  considering  the  simplest 
model  for  nonlinear  disturbance  interactions,  namely  Craik’s  model  [12],  which  provides  a  simplified  account 
of  nonlinear  interactions  via  a  viscous  critical  layer.  Despite  its  shortcomings  from  the  viewpoint  of  systematic 
asymptotics  [14],  this  model  has  been  shown  to  yield  predictions  which  agree  reasonably  well  with  many 
experimental  observations  [20].  In  the  present  context,  it  provides  a  simple  model  to  illustrate  the  combined 
effects  of  nonlinearity  and  randomness  on  the  propagation  of  instability  waves.  Moreover,  by  connecting  the 
amplification  of  subharmonics  to  the  growth  of  phase  coherence,  this  model  suggests  a  link  between  transition 
and  turbulence,  where  the  growth  of  phase  coherence  is  closely  connected  to  energy  transfer.  Fully  developed 
turbulent  flow  exists  because  of  a  balance  between  triad  interactions,  which  promote  phase  coherence,  and 
coupling  to  all  other  modes  which  destroys  it  [23].  The  resonant  triad  model  can  be  understood  from  this 
viewpoint  as  the  first  phase  in  establishing  this  balance. 

We  next  outline  the  extension  of  the  resonant  triad  model  to  more  complex  interactions.  Specifically,  we 
examine  the  interactions  between  more  than  three  instability  waves,  nonlinear  effects  in  a  non-equilibrium 
critical  layer,  and  predictions  for  the  strongly  nonlinear  regime.  These  extensions  raise  the  important  issue 
of  modeling  the  effects  of  higher  order  modes  which  have  not  been  explicitly  resolved  by  the  model  being 
considered.  The  model  which  seems  best  suited  for  engineering  predictions  of  transition  and  could  potentially 
extend  into  the  strongly  nonlinear  regime  is  the  parabolized  stability  equations  (PSE)  [16].  A  preliminary 
discussion  is  given  concerning  how  PSE  can  be  treated  in  a  stochastic  setting. 

Finally,  we  describe  our  proposal  for  a  general,  integrated  turbulence  and  transition  model:  a  generalized 
PSE  restricted  to  a  relatively  small  number  of  modes  with  the  effect  of  the  remaining  modes  accounted  for 
through  modified  damping  and  random  forcing.  The  crux  of  our  proposal  is  that  turbulence  effects  can  also 
be  included  in  the  model  the  same  way.  Such  a  model  would  be  a  compromise  between  transition-sensitized 
Reynolds-averaged  Navier-Stokes  (RANS)  and  large  eddy  simulation  (LES)  or  direct  numerical  simulation 
(DNS). 

2.  Stochastic  aspects  of  the  disturbance  generation  and  evolution  through  the  weakly  non¬ 
linear  stage. 

2.1.  Summary  of  the  receptivity  problem.  The  statistical  analysis  of  transition  in  a  convectively 
unstable  flow  begins  with  the  random  initial  conditions,  which,  barring  uncertainties  or  inaccuracies  in  the 
mean  flow  prediction,  are  the  only  source  of  randomness  in  the  problem.  In  the  Goldstein-Ruban-Zavol’skii 
theory  [15],  [36],  [52],  receptivity  is  analyzed  as  the  excitation  of  a  free  Tollmien-Schlichting  wave  through 
the  interaction  between  a  specified  ambient  unsteady  disturbance  with  surface  roughness  of  known  shape 
and  size.  Provided  that  the  wavenumber  spectrum  of  the  surface  roughness  and  the  characteristic  frequency 
of  the  ambient  disturbance  overlap  respectively  the  wavelength  and  frequency  of  a  free  Tollmien-Schlichting 
wave,  in  low-speed  flows,  this  wave  can  be  excited  and  subsequently  amplified. 

Receptivity  analysis  is  most  often  formulated  as  a  deterministic  problem.  Although  some  stochastic 
aspects  have  been  explored  in  [8],  this  analysis  was  limited  to  evaluating  the  rms  amplitude  of  the  instability 
wave.  By  averaging  over  the  roughness  distribution,  phase  information  is  lost,  and  only  averages  remain. 
This  work  has  been  extended  to  transient  analysis  of  receptivity  by  formulating  the  problem  as  a  stochastic 
differential  equation.  The  result  is  a  description  of  the  growth  of  receptivity  amplitudes  by  means  of  a 
Fokker-Planck  equation  [37].  The  conclusion  from  this  work  for  the  present  purpose  is  simply  this:  because 
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receptivity  is  a  linear  problem,  if  the  phases  in  the  disturbance  environment  are  random,  then  the  phases  in 
the  initial  random  distribution  of  Tollmien-  Schlichting  waves  will  also  be  random. 

While  randomness  of  the  phase  has  no  influence  on  the  linear  evolution  of  the  instability  waves,  there 
exist  cases  in  which  the  disturbance  phase  can  have  a  significant  impact  on  the  nonlinear  evolution  as 
illustrated  later  in  this  section. 

2.2.  Resonant  triad  evolution  in  a  boundary  layer.  The  defects  of  deterministic  transition  pre¬ 
dictions  based  on  linear  amplification  alone  are  well-known:  the  analysis  cannot  account  for  the  rapid  onset 
of  three-dimensional  behavior  as  transition  is  approached,  or  for  the  mean-flow  modification  due  to  the  dis¬ 
turbance  growth.  This  is  because  both  of  these  features  are  due  to  disturbance  nonlinearity  and  the  neglect 
of  nonlinear  effects  on  disturbance  evolution  eventually  invalidates  the  theory.  In  the  following  subsections, 
we  illustrate  the  dependence  of  disturbance  evolution  on  the  relative  initial  phase  of  different  instability 
modes  in  the  simplest  context  to  set  the  stage  for  a  stochastic  analysis  of  nonlinear  disturbance  evolution. 

Many  weakly  nonlinear  analyses  have  been  introduced  in  transition  theory  to  explain  the  onset  of  three- 
dimensional  behavior.  We  consider  the  simplest  theory  of  this  type,  namely  Craik’s  theory  [12]  of  resonant 
interaction  between  a  triad  of  Tollmien-Schlichting  waves  which  shows  how  oblique  subharmonics  can  be 
amplified,  even  if  they  are  linearly  damped. 

2.2.1.  Craik’s  model:  deterministic  triad  interaction  problem.  The  analysis  [12]  begins  with  a 
primary,  two-dimensional  Tollmien-Schlichting  wave  in  a  two-dimensional  parallel  mean  flow 

(2.1)  -  Re{A'Sfo{y)^ax-^} 

which,  together  with  a  pair  of  resonant  subharmonics, 

(2.2)  $1,2  - 

represents  a  special  case  of  a  resonant  triad  of  waves.  Here,  x,y,z  denote  the  streamwise,  normal,  and  the 
spanwise  directions  respectively;  t  denotes  time,  a,/3,  and  u>  denote  the  disturbance  wavenumber  and  fre¬ 
quency  for  a  given  mode,  fa  ( i  —  1, 2, 3)  denote  the  appropriately  normalized  eigenfunctions.  The  disturbance 
amplitudes  At  (i  —  1,2, 3)  are  slowly  varying  functions  of  space, 

1  dA 

(2.3)  Ai  —  Ai(x)  with  |  - - -1  |  <  a 

Ai  ox 

or  of  time, 

1  dA 

(2.4)  Ai  —  Ai(t )  with  |  —  |  <  u> 

depending  on  whether  the  disturbance  evolves  in  space  x  or  time  t.  The  latter  distinction  is  perhaps 
unimportant  for  this  paper.  However,  to  facilitate  comparison  with  previous  deterministic  analyses,  we  will 
treat  only  the  spatial  case  herein. 

Under  the  parallel-flow  assumption,  Craik’s  heuristic  analysis  led  to  the  following  set  of  ordinary  dif¬ 
ferential  equations  for  the  (complex)  amplitudes  of  the  resonant  triad  of  waves  defined  in  Eqs.  (2.1)  and 

(2.2): 

h  =  ah  +  b3b *2 
b,  =  <jb‘>  +  bdb\ 

(2.5)  b3  =  abd  +  e^hb-2 
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where  bi,  (i  —  1,2,3)  denote  normalized  forms  of  the  complex  amplitudes  Ai  and  the  dot  denotes  differen¬ 
tiation  along  the  evolution  coordinate  (i.e.  x  for  the  spatial  disturbance  evolution  considered  here).  In  Eq. 

(2.5) ,  <j  and  a  are  the  normalized  linear  growth  rates  of  the  subharmonic  and  primary  waves  respectively, 
and  the  angle  (j>  models  the  effect  of  detuning  between  the  primary  and  subharmonics.  Eq.  (2.5)  applies 
only  in  the  parallel-flow  limit,  in  which  the  nonlinear  growth  is  assumed  to  occur  on  a  scale  that  is  much 
smaller  than  the  scale  of  the  mean  boundary  layer  growth  (i.e.  nonparallel  scale).  Eq.  (2.5)  also  applies 
more  generally  to  problems  in  which  the  growth  rates  become  weak  functions  of  x  through  small  effects  of 
non-parallism.  This  generalization  will  be  used  in  a  subsequent  case  study. 

An  important  simplification  of  Eqs.  (2.5)  is  provided  by  the  equations  for  parametric  excitation-. 

h  =  abi  +  hb% 
b2  =  &b2  +  b-ibl 

(2.6)  b-i  -  ab3 

wherein  the  nonlinear  term  in  the  equation  for  the  primary  amplitude  b-i  has  been  dropped,  hence  b-s  under¬ 
goes  purely  linear  growth.  This  limiting  case  corresponds  to  the  parametric  region  |  &i  |,  |  f>2  |^|  &3  |  which 
is  typical  of  the  initial  phase  of  nonlinear  growth. 

In  terms  of  amplitude-phase  variables  defined  by 

(2.7)  bi  =  riei$i, 

Eq.  (2.5)  becomes 


r1  —  r2rs  cos  9  +  dri 
f  2  —  rir3  cos  9  +  ar2 
h  —  n r2  cos (4>  -  9)  +  crr-i 

/ci  Q\  a  r'ir-  •  (A.  a\  ri?'3  •  a  r-T'3  •  a 

(2.8)  9  — - sin(®  -9) - sm  9 - sin  9 

r-s  r2  ri 

Observe  that  the  amplitude  evolution  depends,  not  on  the  individual  phases  9i ,  but  only  on  the  single  phase 
variable 

(2.9)  9  =  9i  +  02  -  9j 

Many  special  solutions  and  properties  of  the  triad  equations  can  be  derived  [12],  [48];  however,  for 
the  purposes  of  transition  studies,  a  simple  generic  picture  is  adequate.  In  general,  the  subharmonics  are 
amplified,  eventually  overtaking  the  primary.  During  growth,  the  amplitudes  of  the  subharmonics  become 
equal  (n  =  r2  —  r)  and  a  condition  of  phase-locking  (9  —  0)  develops.  After  the  onset  of  phase-locking,  the 
solution  undergoes  explosive  growth,  leading  to  a  finite-time  singularity  of  the  solution  [48],  Amplification 
can  occur  even  if  the  subharmonics  are  linearly  stable  (cr  <  0),  and  also  even  if  the  primary  is  linearly  stable 
(a  <  0)  [30]. 

This  behavior  is  immediately  apparent  from  the  closed-form  solution  for  the  case  of  parametric  excitation 
Eq.  (2.6): 

h (t)  -  |[6r(0)  +  b2(0 Y]^-1)/^  +  1^(0)  -  b-M*]e~{^~1),,T+&t 

(2.10)  b2(t)  -  |[6i(0)*  +  b2{0)]^ert~1)/a+9t  ~  IMG)*  ~ 
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which  shows  that  the  subharmonics  can  amplify  exponentially  even  if  5  —  0.  Ecp  (2.10)  illustrates  the 
two  features  noted  above:  at  sufficiently  long  time,  amplitudes  of  the  two  subharmonics  become  equal  (i.e. 
I  bi  |  =  |  b‘>  |),  and  bib2  becomes  real  so  that  0  approaches  zero. 

The  parametric  excitation  case  also  reveals  that  amplification  is  phase-dependent:  if  b2(0)  —  —bi(0)  is 
real,  so  that  the  initial  phase  0(0)  =  i r,  then  the  growing  double  exponential  term  does  not  occur  and  we 
have 


h  (t)  =  h  (0)e-(e<T*  _1)/cr+^t 

(2.11)  b2(t)  =  -b1(0)e-<e',-1)/,r+*t. 

In  this  case,  the  subharmonics  will  actually  decay  in  the  absence  of  linear  growth. 

We  note  that  the  dependence  of  subharmonic  amplification  on  its  phase  relative  to  the  primary  wave 
is  well-known  in  the  context  of  mixing-layer  instabilities  [30].  Analogous  behavior  in  wall-bounded  flows, 
albeit  not  equally  emphasized  in  the  literature,  has  been  measured  [39]  and  theoretically  predicted  [53]. 
Unfortunately  however,  the  effects  of  initial  phase  lag  on  transition  location  have  not  been  explored  in  detail 
as  yet.  A  preliminary  attempt  using  admittedly  crude  modeling  is  made  below  (see  Fig.  2.8).  Further 
examination  using  both  nonlinear  PSE  and  wind  tunnel  experiments  is  planned  for  the  near  future. 

The  general  features  of  triad  evolution  are  depicted  in  Fig.  2.1  which  were  obtained  from  numerical 
integration  of  Eqs.  (2.5)  using  a  fourth-order  Runge-Kutta  method  starting  from  initial  values  of  the  phase 
0(0)  =  iV7r/4,  N  —  1,  •  •  •  16.  In  all  cases,  initial  amplitudes  were  held  fixed  at  n(0)  =  r*2(0)  =  0.1,  rg(0)  =  1.0 
with  initial  growth  rates  a  —  0.  (i.e.  linearly  neutral  subharmonies),  and  a  —  0.1.  Figure  2.1  shows 
both  the  evolution  of  the  phase,  in  which  “locking”  is  clearly  evident,  and  the  evolution  of  the  primary 
and  subharmonic  amplitudes;  a  strong  connection  between  subharmonic  amplification  and  phase  locking  is 
suggested  by  a  comparison  between  the  phase  and  amplitude  evolution. 


Fig.  2.1.  Effect  of  initial  phase  on  subharmonic  amplification  in  the  nonlinear  evolution  of  a  resonant  triad.  The 
downstream  development  of  the  phase  variable  9  and  of  the  primary  and  subharmonic  amplitudes  are  shown  for  variable  initial 
9.  Parameters  in  Eq.  (2.5)  are:  a  —  0.1  and  5  —  0;  the  dephasing  factor  —  0.  Initial  conditions  are:  primary  amplitude 
7*3(0)  =  1.0  and  subharmonic  amplitudes  7*i(0)  =  7*2(0)  =  0.1.  The  initial  phase  0(0)  is  uniformly  distributed  from  w  to  w. 
The  governing  equations  are  fully  nonlinear  triad  interactions ,  Eq.  (2.5). 
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Observe  that  even  though  the  triad  system  Eq.  (2.5)  supports  an  equilibrium  solution  of  the  form  9  —  n, 
that  solution  is  unstable,  with  the  only  long-time  attractor  being  the  phase-locked  solution  9  —  0  as  observed 
earlier  by  [14]  in  a  similar  context  using  high  Reynolds  number  asymptotic  theory. 

Fig.  2.2  shows  the  same  calculation  for  the  linear  theory  of  parametric  excitation.  The  picture  is 
naturally  considerably  simpler  since  the  primary  only  undergoes  linear  growth.  The  curve  with  9(0)  —  tt 
maintains  its  initial  value  througout.  Therefore,  unlike  the  fully  coupled  case  illustrated  by  Fig.  2.1,  the 
subharmonic  amplitude  in  this  case  continues  to  decay  during  its  evolution. 


Fig  .  2.2.  Effect  of  initial  phase  on  subharmonic  amplification  in  the  case  of  parametric  excitation.  Disturbance  parameters 
and  initial  conditions  are  as  in  Fig.  2.1. 

2.2.2.  Stochastic  formulation  of  triad  evolution:  probability  density  function  approach. 

The  analysis  of  subharmonic  growth  suggests  a  role  for  a  statistical  theory  of  resonant  interactions.  Since 
the  growth  of  subharmonic  disturbances  depends  on  their  initial  phase,  an  exact  description  will  require  a 
knowledge  of  the  initial  phase;  but  if  these  initial  conditions  are  referred  to  a  receptivity  calculation,  then 
as  Sect.  2.1  suggests,  even  if  the  initial  disturbance  amplitudes  were  known  with  relative  precision,  the 
initial  phase  distribution  could  not  be  established  in  general.  Thus,  the  initial  conditions  for  the  theory  of 
subharmonic  growth  are  best  described  by  a  probability  density  function  (PDF)  for  initial  amplitudes  and 
phases.  From  this  viewpoint,  the  best  description  of  subharmonic  growth  is  given  by  the  evolution  of  this 
PDF. 

For  the  sake  of  brevity,  we  formulate  the  evolution  equation  for  the  case  of  equal  subharmonic  amplitudes, 
ri  =  r-s  —  r.  This  condition  is  preserved  by  the  equations  of  motion,  and  as  noted  previously,  the  system 
eventually  evolves  to  this  condition  in  any  case.  Applying  standard  methods  [31],  one  finds  that  the  joint 
probability  density  P(r,  r-i,9;  t)  of  amplitudes  r,  and  the  phase  variable  9  satisfies  the  first  order  partial 
differential  equation 

(2T2)  f +  ±{FrP)  +  ±.{F9P)  +  ±iF.P)  =  0 

where 


Fr  —  rr-j  cos  9  +  dr 
F'z  —  r 2  cos {(/>  —  9)  +  0T3 
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(2.13) 


Fo  —  —  sin ((f)  —  9)  —  2 r3  sin  9. 

r-3 

Closed  equations  for  the  moments  of  P  cannot  be  obtained,  since  moments  of  any  given  order  are  found  to 
depend  on  moments  of  higher  order. 

The  initial  PDF  is  taken  to  be  independent  Gaussian  distributions  for  amplitudes,  and  initially  uniformly 
distributed  phase, 

(2.14)  P(r,  r3,9;  0)  =  J-N(r  \  r0,a0)N(r3  \  r3t0,  a3) 

Ztt 

where 

(2.15)  N{x  |  xo,  a)  =  1  e-(»~«o )2/2^. 

V  27 r<7“ 


Since  Eq.  (2.12)  is  a  linear  first-order  partial  differential  equation  in  a  conservative  form,  numerical 
integration  is  straightforward.  The  method  used  for  solving  this  equation  numerically  is  described  in  detail 
in  Appendix  I.  The  numerical  integration  reproduced  the  qualitative  features  of  the  deterministic  analysis  as 
described  earlier  via  Figs.  2.1  and  2.2.  The  marginal  probability  densities  for  subharmonic  amplitude  and 
phase, 

poo  p2ir 

P{r3,t)=  dr  d9  P(r,r3,9;t) 

Jo  Jo 

poo  poo 

(2.16)  P(9,t)=  dr  dr3  P{r,r3,9-,t) 

J  o  Jo 

are  shown  for  the  case  of  parametric  excitation  in  Fig.  2.3. 

This  case  is  chosen  for  preliminary  code  verification  because  the  primary  amplitude  probability  density 
must  be  constant  as  a  function  of  the  centered  normalized  variable  7  -  (r  —  (r))/cr.  The  primary  amplitude 
PDF  remains  Gaussian  and  unchanged  after  this  normalization.  Fig.  2.3  shows  that  the  subharmonic 
amplitude  becomes  bimodal  as  time  increases.  The  bimodality  arises  because  paths  exist  for  which  the 
subharmonic  amplitude  decays  at  first.  This  M-shaped  PDF  can  be  compared  with  the  measurements  of 
[29]  in  the  analogous  case  of  a  transitioning  mixing  layer.  Fig.  2.3  also  presents  another  picture  of  the  strong 
tendency  to  phase  locking. 

For  quantitative  validation  of  the  method,  we  compare  the  results  with  some  partially  analytic  results 
available  for  the  case  of  parametric  excitation.  A  straighforward  calculation  shows  that  the  distribution  of 
phase  is  related  to  the  initially  uniform  distribution  ip  through  the  equation 


(2.17) 


d$  -  -2  R — Tj- 
sin' 


1  —  cos  ip 

ip  +  (cos  ip  —  1)2R2 


dip 


where  ip  is  uniformly  distributed,  and  the  factor  R  is  the  double  exponential  of  Eq.  (2.10),  i.e., 


(2.18)  R  —  exp[—  exp(<rt  —  1)  +  at]. 

The  marginal  probability  density  function  of  the  phase  is  compared  with  the  exact  result  in  Fig.  2.4  .  The 
agreement  is  satisfactory.  Analytical  expressions  for  the  amplitude  PDF’s  cannot  be  derived;  however,  the 
mean  and  variance  of  the  primary  and  subharmonic  are  easily  evaluated  numerically  in  this  case.  The  results 
are  compared  with  the  predictions  of  the  PDF  code  in  Fig.  2.4.  The  agreement  is  satisfactory  except  that 
the  predictions  of  the  subharinonic  standard  deviation  become  inaccurate  near  the  end  of  the  simulation. 
This  inaccuracy  may  reflect  the  low  order  accuracy  of  the  finite  volume  algorithm  used  in  the  code. 
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Fig.  2,3.  Subharmonic  amplitude  and  phase  PDF’s  for  different  Re  under  parametric  excitation.  Initial  phase  is  uniformly 
distributed  from  tv  to  n;  initial  primary  amplitude  is  Gaussian  with  mean  1  and  standard  deviation  0.1,  and  equal  initial 
subharmonic  amplitudes  are  Gaussian  with  mean  0.1  and  standard  deviation  0.01.  Disturbance  parameters  are  as  in 

Fig.  2.1.  hs  Re  increases,  the  initially  Gaussian  subharmonic  amplitude  PDF  becomes  bimodal  and  the  initially  uniform  phase 
PDF  becomes  concentrated  at  0  —  0. 


Fig.  2.4.  Validation  of  numerical  solution  of  joint  PDF  evolution  equation  for  parametric  excitation  case:  initial  conditions 
are  as  in  Fig.  2.3.  Left  -  comparison  between  exact  (dots)  phase  evolution  and  numerical  solution  (lines)  of  joint  PDF  equation. 
Right  -  comparison  between  exact  (symbols)  and  computed  (lines)  amplitude  statistics.  Disturbance  parameters  o,a,4>  are  as 
in  Fig.  2.1. 

2.3.  Stochastic  form  of  triad  evolution:  particle  method.  Although  direct  solution  of  the  PDF 

evolution  equation  Eq.  (2.12)  is  feasible,  the  computation  time  required  is  excessive.  Even  with  the  as¬ 
sumption  that  the  amplitudes  of  the  two  subharmonic  modes  are  equal,  integration  up  to  the  time  that  the 
subharmonics  exceed  the  primary  required  several  days  on  a  workstation.  As  noted  in  Ref.  [31],  the  main 
difficulty  in  integrating  any  PDF  evolution  equation  is  dimensionality:  updating  a  PDF  with  d  variables 
with  N  degrees  of  freedom  each  requires  about  Nd  operations  per  time  step. 

The  particle  method  of  Pope  [31]  mitigates  this  difficulty  at  the  expense  of  a  coarser  description  of 
the  PDF.  In  the  present  problem,  the  particle  method  can  be  formulated  as  follows.  Consider  the  general 
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problem  with  independent  amplitudes  rx,r2,r3.  Replace  the  initial  PDF  by  a  sum  of  delta  functions, 


P(ri,r2,r3,9\0)  «  ^  P(ni,ri2,n3,ne)  x 

m,n2,n3,n^ 

(2.19)  <S(n  -  rm)5(r2  -  rn2)5(r3  -  rn3)J(0  -  0n„) 

where  the  points  rrii,9ris  at  which  the  PDF  is  sampled  are  arbitrary,  but  deterministic.  These  points  could 
be  chosen  based  on  a  Gaussian  quadrature  scheme,  for  example.  Note  that  if  the  variables  ry,  9  are  assumed 
to  be  independent  initially,  then 

(2.20)  P{ni,n2,n3,ne)  —  P{n1)P(n2)P(n3)P(9) 

Next,  the  deterministic  evolution  equations  are  integrated  for  an  ensemble  of  initial  conditions  n(0)  — 
rni,r2(0)  —  rn2,r3(0)  —  rns,9(Q)  —  9ne.  The  PDF  at  time  t  is  approximated  by  the  discrete  sum 

P(ri,r2,r3,9-,t)  w  ^  P{nun2,n3,n9)  x 

m,n2,n3,n9 

(2.21)  d(n  -  rni(t))S(r2  -  rn.Jt))6(r3  -  rri3(t))6(9  -  9n„(t)). 

This  scheme  is  not  a  Monte  Carlo  method,  because  the  initial  PDF  is  not  sampled  randomly;  instead,  the 
sampling  is  deterministic  but  weighted  according  to  the  specified  initial  PDF.  An  important  feature  of  the 
method,  which  is  evident  from  Eq.  (2.21),  is  that  it  conserves  total  probability. 

To  assess  the  accuracy  of  this  particle  method,  the  statistics  shown  in  Fig.  2.4  were  recomputed  and 
compared  with  the  exact  results  in  Fig.  2.5.  Evidently,  the  agreement  is  much  more  satisfactory. 


X 


Fig.  2.5.  Comparison  between  theoretical  and  computed  amplitude  statistics  for  parametric  excitation:  tines  -  computation 
(particle  method),  symbols  -  theory.  Initial  conditions  and  disturbance  parameters  0,  o,<f>  are  as  in  Fig.  2.1. 

One  must  note,  however,  that  one  of  the  advantages  of  the  probability  density  function  approach  is  that 
it  is  robust  with  respect  to  any  finite-time  singularities  observed  in  individual  solutions  of  the  deterministic 
evolution  equations.  The  particle  method,  on  the  other  hand,  requires  one  to  remove  any  trajectories  that 
approach  a  singularity  (and  to  account  for  this  removal  in  some  reasonable  manner).  However,  this  is  not 
necessarily  a  significant  limitation  as  the  singularities  indicate  limitations  of  the  underlying  (deterministic) 
model  and,  hence,  are  either  delayed  or  bypassed  altogether  in  a  real  application  after  the  model  is  suitably 
generalized  to  include  the  relevant  physics  neglected  previously. 
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2.4.  Stochastic  triad  evolution:  case  study  and  comparison  with  experiment.  The  particle 
method  was  applied  to  examine  stochastic  aspects  of  the  triad  evolution  in  a  flat-plate  boundary  layer.  In 
particular,  the  case  of  a  symmetric  resonant  triad  which  had  previously  been  examined  in  Ref.  [53]  in  a 
deterministic  setting  was  considered.  This  triad  was  comprised  of  a  two-dimensional  fundamental  wave  with 
a  dimensionless  frequency  F  —  115  x  10-6  and  its  subharmonics  with  a  spanwise  wavenumber  given  by 
j3/Re  —  0.22  x  10-6.  (For  a  definition  of  these  dimensionless  parameters,  the  reader  is  refered  to  [53].) 

To  allow  for  the  effects  of  the  growth  of  the  mean  boundary  layer,  as  well  as  the  effects  of  small  detuning, 
the  triad  equation  set  Eq.  (2.5)  was  replaced  by  the  slightly  more  general  system 

bi  -  a{x)bi  +  5i,23(a0&2&3 
b2  -  a(x)b2  +  52,i3(^)^^3 

(2.22)  b3  -  cr(x)b3  +  S3, i2{x)hb2 

where  the  x  dependence  of  the  linear  growth  rates  a  and  <7  and  the  nonlinear  interaction  coefficients  reflect 
the  effects  of  weak  mean-flow  non-parallelism.  The  interaction  coefficients  S[jh  can  be  complex  to  allow 
for  detuning  effects.  The  evolution  variable  is  streamwise  distance,  measured  in  this  case  as  a  downstream 
Reynolds  number.  This  system  can  also  be  written  in  terms  of  amplitude  and  phase  variables  as 

h  -  oti  +  Ri,23(x)r2r-i  cos (9  +  <f>i(x)) 
f2  -  ar2  +  R2,3i{x)r2ri  cos (9  +  4>2{x)) 
h  —  crr-i  +  R3,i2(x)rir2  cos (03(a;)  -  9) 

(2.23)  9  =  R3,i2—  sin (03 (x)  -9)-  Jfc,3i  —  sin(0  +  4>2{x)) 

T  3  r-2 

where  we  have  set 

(2-24)  SiJk  =  Rijke**' . 

The  values  of  the  nonlinear  coefficients  Sijk  as  functions  of  the  Reynolds  number  were  based  on  Table 
1  of  [53].  The  initial  PDF  is  given  in  Eq.  (2.14). 

The  results  are  shown  in  Fig.  2.6.  The  right-hand  graph  can  be  compared  with  Fig.  2  of  [53],  The 
subharmonic  amplitudes  are  initially  unequal;  however,  the  tendency  of  these  amplitudes  to  equalize  is 
evident.  The  subharmonics  overtake  the  primary  at  a  Reynolds  number  of  about  700.,  in  agreement  with  the 
analysis  of  [53].  After  this  point,  the  entire  system  undergoes  explosive  growth  leading  to  a  singularity  in  all 
three  amplitudes.  The  left-hand  graph  shows  the  corresponding  phase  evolution.  Both  the  close  connection 
between  the  development  of  phase  locking  and  the  onset  of  explosive  nonlinear  growth  and  the  effect  of  initial 
phase  on  subharmonic  growth  which  were  observed  in  the  previous  model  problems  also  exist  in  this  more 
realistic  example. 

In  order  to  illustrate  the  connection  between  subharmonic  amplification  and  transition,  the  conditions 
analyzed  in  [53]  were  modified  slightly.  In  Fig.  2.7,  the  initial  amplitudes  of  the  two  subharmonics  are 
one-half  the  initial  primary  amplitude.  Initial  standard  deviations  of  the  amplitudes  are  set  at  ten  percent 
of  the  mean.  The  effect  of  initial  phase  difference  on  subharmonic  growth  is  considerably  enhanced  by  this 
change  of  initial  amplitude. 

To  connect  this  analysis  to  transition,  Fig.  2.8  shows  the  probability  that  the  subharmonic  amplitude 
exceeds  the  local  primary  amplitude.  Although  the  crossover  location  is  not  the  same  as  the  transition 
onset  location  (which  cannot  be  predicted  within  the  limited  framework  of  resonant  triad  interaction),  it 


-  #1,23  sin (6*  +  0i  (x)) 
r  1 
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Fig.  2.6.  Phase  locking  and  subharmonic  amplitude  equalization  during  triad  evolution  in  the  problem  analyzed  in  Ref. 
[53].  Left  -  phase  evolution  as  a  function  of  Reynolds  number.  Right  -  amplitudes  of  primary  and  subharmonics  as  functions 
of  Reynolds  number.  The  right-hand  figure  can  be  compared  with  Fig.  2  of  Ref.  [53]. 


Fig.  2.7.  Effect  of  initial  phase  on  triad  evolution.  Primary  and  subharmonic  amplitudes  are  shown  as  functions  of 
local  Reynolds  number.  All  triad  evolution  parameters  are  based  on  Fig.  2  of  Ref.  [53],  except  that  the  initial  subharmonic 
amplitudes  are  both  equal  to  0.5  times  the  initial  primary  amplitude. 


does  harbinger  the  potential  onset  of  stronger  nonlinear  interactions  and,  hence,  the  approach  of  transition 
(assuming,  of  course,  that  the  primary  amplitude  involved  is  significantly  large  as  well). 

With  a  suitable  generalization  of  the  underlying  model,  then,  a  stochastic  transition  theory  will  indicate 
the  probabilistic  spread  in  transition  location  as  a  function  of  uncertainties  in  the  initial  disturbance  char¬ 
acteristics.  In  that  spirit,  Fig.  2.8  may  be  viewed  as  a  measure  of  variation  in  transition  onset  probability 
with  downstream  location. 

A  qualitative  comparison  with  experimental  results  [39]  is  shown  in  Fig.  2.9.  Here,  the  subharmonic 
amplitude  at  selected  downstream  locations  is  plotted  as  a  function  of  initial  phase  difference  between  the 
primary  and  subharmonic  waves. 
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Fig.  2.8.  Transition  onset  probability  based  on  subharmonic  growth  as  function  of  local  Reynolds  number.  The  conditions 
analyzed  are  the  same  as  in  Fig.  2. 7.  Transition  onset  is  assumed  to  occur  when  the  primary  and  subharmonic  amplitudes  are 
equal.  The  spread  in  transition  onset  location  is  due  to  uncertainty  of  the  initial  phase  and  amplitude. 


Fig.  2.9.  Subharmonic  amplitude  as  a  function  of  initial  phase  difference  between  subharmonic  and  primary  modes:  left 
-  theory ,  right  -  experiment  of  [39].  The  triad  parameters  are  somewhat  different,  therefore  the  comparison  is  qualitative. 


3.  Extension  of  stochastic  formulation  to  more  general  nonlinear  disturbance  evolution 
models.  By  accounting  for  the  rapid  onset  of  three-dimensionality  in  transitional  flows,  the  resonant  triad 
model  overcomes  one  of  the  important  limitations  of  linear  stability  theory;  this  was  actually  one  of  the 
motivations  behind  the  original  development  of  the  theory.  Nevertheless,  the  resonant  triad  model  cannot 
be  used  to  predict  the  onset  of  transition  in  general  flows.  One  obvious  problem  is  that  by  limiting  the 
analysis  to  a  system  of  only  three  interacting  modes,  resonant  triad  theory  can  describe  neither  a  wide-band 
spectrum  of  initial  disturbances  nor  the  generation  of  a  wide-band  spectrum  through  nonlinear  interactions  as 
transition  proceeds.  Moreover,  like  linear  theory,  it  predicts  an  indefinite  growth  of  the  disturbances,  which 
ultimately  contradicts  the  starting  assumption  of  weak  nonlinearity.  Finally,  transition  may  not  necessarily 
occur  as  a  consequence  of  resonant  triad  interaction:  there  exist  other,  qualitatively  different  scenarios  such 
as  oblique  or  fundamental  breakdowns  and  wave- vortex  interactions,  which  may  initiate  transition  and  which 
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may  require  a  different  type  of  modeling. 

We  will  therefore  discuss  other  deterministic  disturbance  evolution  models  and  their  implications  for 
statistical  modeling. 

3.1.  Overview  of  deterministic  models.  The  simplest  way  to  improve  the  resonant  triad  model 
is  to  increase  the  number  of  modes  in  the  analysis.  Although  this  type  of  model  remains  restricted  to 
weak  nonlinearity,  the  growing  disturbances  can  be  spread  over  a  larger  number  of  modes;  the  consequent 
possibility  that  the  amplitude  of  any  one  mode  can  saturate  might  extend  the  limits  of  applicability  of  the 
assumption  of  weak  nonlinearity.  Evidence  supporting  this  conjecture  comes  from  the  analysis  of  a  five- wave 
resonant  system  in  [53],  in  which  it  is  found  that  the  effect  of  additional  modes  is  to  delay  the  onset  of 
explosive  nonlinear  growth. 

Standard  perturbation  methods  can  be  used  to  derive  governing  equations  for  more  complex  interactions 
of  wave  amplitudes.  The  simplest  generalization  of  the  resonant  triad  model  is  [53] 

(3.1)  (-7^  +  W^i—  +  At  -  7 (A(  +  y^hu+kSijkAjAk 

'  j,k 

where  the  indices  l,k,j  vary  over  a  finite  set  of  modes.  In  addition  to  describing  the  interactions  of  a  larger 
set  of  modes,  Eq.  (3.1)  allows  for  off-resonant  interactions  through  the  variation  of  the  phase  factor  hij+k- 
Under  strongly  off-resonant  conditions,  the  various  modes  can  still  communicate  through  the  mean-flow 
correction  (A0)  which  is  retained  in  this  model.  Of  course,  the  significant  distortion  of  mode  shapes  (and 
concomitant  modification  of  the  growth  rates  7 ;)  under  a  strong  enough  mean-flow  correction  is  not  explicitly 
accounted  for  in  the  model.  The  model  also  allows  spanwise  modulations  of  the  modal  amplitudes. 

In  the  framework  of  the  same  perturbation  scheme,  the  equation  which  incorporates  the  next  order 
corrections  is  [53] 


f  d  t  w  d  w  d  if  0'2u>i  d 2  i  d2u>i  d2  t  d2u>i  d2  \  ]  A 

+  Wl'ldfc+  2,l"dz  ~  ~  2  V  da2  dx~  +  dad/3  dxdz  +  d/P  dz 2  )  }  1 

—  "^2,  ^l-jkAjAk  +  (Cljl-fa  +  ^  hij+k  +  T:  CljkrhlJ+k+rAjAkA , 


This  equation  exhibits  higher  order  nonlinearities  through  the  cubic  coupling  terms,  a  more  complex, 
amplitude-dependent  dispersion  through  the  terms  proportional  to  Cyl  and  Cj?l,  and  dependence  on  the 
spatial  derivatives  of  the  modal  amplitudes.  The  definitions  of  all  the  coefficients  are  given  in  the  original 
reference  [53], 

The  question  arises  whether  the  higher  order  models  like  Eq.  (3.2)  represent  an  essential  improvement 
over  the  lowest  order  approximation  Eq.  (3.1).  The  answer  is  that  this  entire  class  of  nonlinear  theories 
has  fundamental  limitations  which  are  not  addressed  by  formulating  higher  order  theories.  As  mentioned  in 
the  Introduction,  this  entire  class  of  models  can  be  criticized  from  the  viewpoint  of  more  systematic  high 
Reynolds  number  asymptotics,  because  these  theories  do  not  properly  reduce  to  linear  theory  in  the  far 
upstream  limit.  In  fact,  as  noted  by  [53],  the  ‘self-consistent  assignment  of  harmonics  at  x  —  0’  constitutes 
a  difficulty  for  the  theory.  Moreover,  evaluation  of  the  coefficients  in  terms  of  linear  eigenfunctions  shows 
that  their  values  are  dominated  by  contibutions  from  nonlinear  critical  layers.  This  observation  suggests  an 
important  role  for  the  nonlinear  development  of  the  critical  layers. 

Both  of  the  above  considerations  have  led  to  a  re-evaluation  of  triad  and  generalized  triad  theories  from 
the  standpoint  of  systematic  high  Reynolds  number  asympotic  expansions.  A  representative  result  is  the 
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integro-differential  equation  of  Goldstein  and  Lee  [14] 


dA  4  2 . 

~dx~  5A+5* 


(x  —  xi)2A3(xi)A*(2x  —  xi)dxi 


dA2 

dx 


/x  rxi 

/  K\{x  |  xi, x2)A(xi)A(x2)A*(xi  +  X2  —  2x)dx2dxi 

-OO  J  —  OO 

/x  rx  i 

/  [2(x  -  Xi)3A3(xi)A(x2)A*(2xi  +  x2  -  2x)  + 

-OO  v  — OO 


(3.3) 
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K2{x  |  xx,x2)A(xi)A3(x2)A*(xi  +  2x2  -  2x)]dx2dxi 

a;  /*#!  /*a72 

+i  I  /  /  iv3(x  |  xi,x2,x3)A(xi)A(x2)A(X3)A*(xi  +  x2  +  x3  -  2x)dx3dx2dxi. 
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In  this  equation,  the  primary  amplitude  is  A3,  and  the  subharmonic  amplitudes  are  Ai  —  A2  —  A.  The 
coefficients  Kj  and  kernel  functions  Ki  are  known  in  terms  of  closed-form  expressions.  All  of  the  relevant 
definitions  are  given  in  [14]. 

The  form  of  the  integral  terms  on  the  right  sides  of  Eq.  (3.3)  show  that  the  solution  matches  the  linear 
solution  in  the  limit  x  — >  —  oo.  Achieving  the  correct  matching  was  indeed  the  goal  of  the  entire  analysis. 
It  is  reassuring  that,  in  spite  of  the  differences  in  the  form  of  the  nonlinear  term  involved,  solutions  of  Eq. 
(3.3)  in  fact  are  qualitatively  analogous  to  those  of  Craik’s  model.  There  is  an  initial  phase  of  parametric 
amplification  of  the  subharmonics  followed  by  explosive  nonlinear  growth  of  all  three  modes.  The  same 
sensitivity  to  initial  relative  phase  also  exists;  the  nonlinear  growth  phase  of  the  evolution  can  be  delayed 
for  certain  values  of  initial  phase  [14]. 

A  simpler  equation  than  Eq.  (3.3)  has  been  derived  in  [28]: 


(3.4) 


^Aq  _ 

— J —  —  KoAq 

dx 

dA  37 tR3  .  , .  \ 

-j-  —  kqi, A  +  —  — —  iA  Ao  +  iMA{x) 
dx  XU  A 


A(x  —  x')  |2 


The  linear  growth  of  the  primary  suggests  purely  parametric  excitation,  but  the  subharmonic  evolution 
equation  contains  both  an  algebraic  nonlinearity  similar  to  the  Craik  model  (attributed  to  a  viscous  critical 
layer),  and  a  nonlinear  self- interaction  term.  Indeed,  after  dropping  the  self-interaction  term,  Eq.  (3.4) 
reduces  to  a  special  case  of  Eq.  (2.6).  Again,  all  of  the  notation  in  Eq.  (3.4)  is  defined  in  [28].  Analogous 
theories  for  other  types  of  wave  interactions  such  as  wave-vortex  interaction,  oblique  mode  breakdown,  and 
phase-locked  interactions  have  also  been  developed;  see,  for  instance  [10],  [50]. 

While  the  asymptotic  theories  have  aided  greatly  in  our  understanding  of  the  various  transition  scenarios, 
each  of  them  pertains  to  some  specific  kind  of  dominant  physical  balance.  Also,  in  some  cases,  the  underlying 
asymptotic  behavior  is  only  established  at  unpractically  large  Reynolds  numbers.  Hence,  such  models  are 
perhaps  not  well-suited  for  engineering  prediction  of  transition.  Another  transition  model,  which  accounts 
for  mean-flow  nonparallelism  and  disturbance  nonlinearity  in  a  composite  sense,  is  the  PSE  model  introduced 
by  Herbert  [4],  [16].  Albeit  not  as  rigorous  as  the  high  Reynolds  number  asymptotic  theories,  it  is  better 
suited  as  a  general  purpose  transition  prediction  tool  and  has  been  shown  to  yield  surprisingly  accurate 
predictions  through  somewhat  past  the  rise  of  the  mean  skin-friction  curve. 

To  formulate  the  PSE  for  a  two-dimensional  incompressible  mean  flow  Uq(x,t/),  the  disturbance  field 
u(x,  y,  z,  t)  representing  both  velocity  and  pressure  is  written  as  a  Fourier  series 


(3.5) 


u  =  ^Ara,n(x,7/)ei^™t) 

m,n 
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where  the  complex  amplitudes  A„tjn  are  further  decomposed  into  a  slowly  varying  shape  function  and  a 
rapidly  oscillating  phase 


(3.6) 


A-m,n(x,y)  -  '&m,n{x,y)exp 


i  I  (y  ^)dx 

.  J  Xq 


The  wavenumbers  aTOjri.  are  determined  so  that  the  downstream  evolution  of  the  complex  amplitudes  WTOjn  is 
as  slow  and  smooth  as  possible.  While  there  are  many  alternatives  to  implement  this  constraint,  a  common 
one  is  the  integral 

f  d** 

(1.7)  =  0 


where  the  asterisk  denotes  complex  conjugation.  The  above  constraint  both  renders  the  decomposition  in  Eq. 
(3.6)  unique  and  ensures  the  approximate  validity  of  the  marching  approximation  to  the  originally  elliptic 
Navier-Stokes  equations.  The  details  of  this  procedure  are  discussed,  for  example,  in  [6]. 

The  final  PSE  equations  can  be  expressed  in  the  general  form 

/o  o\  n  _i_  A  dtt'in.n  _  d~'i!m,rl 

V01*0/  L/ra,n^ra,n  t  /'ra,n  i  L>m,n  —  v  m,n  ^^2  ' 

Here  the  matrices  ATO>r{,  •  •  •  Dm>r(  depend  on  the  complex  wavenumber  am,n  and  its  streamwise  derivatives. 
The  nonlinear  forcing  term  fm>H  can  be  written  as  a  sum  of  multilinear  forms: 


(3.9) 


fm,n  =  W(*>V¥) 


i+k=m,j+l=n 


exp 


i  I  &i,j  {x  )  "  oc /.■  i  (  j.  )  kXjfi  j { ( ,r  )  (lx 
.  J  Xo 


In  what  follows,  the  indices  on  the  multilinear  forms  B  and  its  arguments  will  be  dropped  to  simplify  the 
notation. 

The  mode  ^0,0  represents  the  mean  flow  perturbation.  The  evolution  equation  for  ’J'o.o 


(3.10) 


Do, 0^0,0 


d^0,o  ,  R  c^^o.o 

W),o — W - r  tS0iQ 


dx 


dy 


52^o,o  f 

—  v0,0 — 7y, - r  Io,0 

dyz 


contains  the  mean  Reynolds  stress  gradients  fy_y.  Hence,  the  equation  for  the  total  mean  velocity  field 
U  =  Uq  +  ^0,0  (where  the  sum  is  restricted  to  velocity  components  only),  is  governed  by 


(3.11) 

where 


d2\J 

U-VU--VP  +  F— 

dy 


+  fo,o 


(3.12)  fo,o=  B[*TOl„,V*_mi_„]. 

The  PSE  equations  provide  a  closure  for  the  term  fy  0  since  evolution  equations  for  are  given  by  Eq. 

(3.8). 

The  PSE  equation  Eq.  (3.8)  is  a  type  of  triple  decomposition  in  which  the  fluctuations  are  projected  onto 
components  which  have  known  evolution  equations  rather  than  being  separated  into  steady  and  unsteady 
parts  as  in  Reynolds  averaging.  The  terms  fTOi„  with  ( m,n )  ^  (0,0)  which  also  arise  in  the  PSE  can  be 
understood  as  generalized  Reynolds  stress  gradients. 
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PSE  calculations  in  this  form  have  typically  been  carried  out  only  for  a  short  distance  past  the  point  at 
which  the  skin  friction  begins  to  rise.  Eventually,  the  iterative  numerical  process  (which  includes  nonlinear 
iterations  coupled  within  the  interations  required  to  determine  the  functions  aTOi„  in  Eq.  (3.6))  fails  to 
converge  before  laminar  breakdown  is  completed.  To  some  extent,  the  nonconvergence  may  be  delayed  by 
treating  the  nonlinear  ‘source7  term  in  an  implicit  fashion  ([2]).  However,  the  validity  (or  at  least  accuracy) 
of  the  PSE  model  itself  becomes  questionable  in  the  laminar  breakdown  region.  This  region  involves  a 
rapid  adjustment  of  the  mean  flow  from  laminar  to  turbulent  flow,  typically  over  a  rather  small  number  of 
wavelengths.  This  means  that  the  terms  in  d2 /dx2  are  not  necessarily  negligible,  and  that  the  decomposition 
of  the  fluctuations  into  “wave”  and  “amplitude”  components  becomes  suspect.  Only  extensive  testing  will 
show  whether  or  not  the  PSE  could  be  applied  into  the  breakdown  region.  Therefore,  from  the  standpoint  of 
integrated  transition  and  turbulence  modeling,  it  will  be  desirable  to  switch  over  to  a  conventional  turbulence 
model  somewhere  in  the  breakdown  region. 

3.2.  Stochastic  formulation:  PDF  evolution  equations  for  modal  amplitudes.  We  now  con¬ 
sider  how  to  recast  the  deterministic  models  from  Sect.  3.1  in  stochastic  form,  specifically  as  a  closed 
evolution  equation  for  the  joint  probability  density  function  of  modal  amplitudes.  These  models  fall  into  a 
natural  hierarchy,  which  progresses  from  pui'ely  algebraic  nonlinearity  (Eq.  (3.1)),  to  integral  nonlinearity 
(Eq.  (3.3)),  to  a  nonlinear  partial  differential  equation  (the  PSE  system  Eq.  (3.8)).  As  already  noted  at 
the  beginning  of  this  section,  except  for  requiring  the  resolution  of  a  PDF  with  a  very  large  number  of 
independent  modal  amplitudes,  there  is  no  difficulty  in  formulating  a  closed  PDF  evolution  equation  cor¬ 
responding  to  systems  like  Eq.  (3.1)  with  purely  algebraic  nonlinearity.  However,  the  formulation  of  PDF 
evolution  equations  for  systems  with  noil-algebraic  nonlinearities  like  Eq.  (3.3)  and  Eq.  (3.8)  will  present 
closure  difficulties. 

An  expanded  set  of  modes,  as  suggested  in  Eq.  (3.1),  is  easily  accomodated  in  the  stochastic  formulation, 
at  least  in  principle.  The  only  problem  with  an  expanded  set  of  modes  is  the  practical  problem  of  solving 
a  PDF  evolution  equation  containing  a  large  number  of  modes.  However,  deriving  evolution  equations  for 
the  joint  probability  density  function  of  modal  amplitudes  corresponding  to  either  high-Reynolds  number 
asymptotic  theories  or  the  PSE  model  introduces  a  new  difficulty  because  these  equations  all  couple  infor¬ 
mation  from  different  points,  either  through  the  spatial  derivatives  which  occur  in  Eqs.  (3.1)  and  (3.2),  or 
through  the  spatial  integral  in  Eq.  (3.3).  This  coupling  between  different  points  makes  it  impossible  to 
derive  an  equation  for  the  single-point  probability  density  without  closure  assumptions  as  described  below. 

In  general  terms,  given  an  evolution  equation  for  a  random  field  A(x) 

(3.13)  A  —  /(A,  VA) 

where  the  dot  denotes  differentiation  with  respect  to  the  evolution  variable,  the  probability  density  V{A) 
satisfies  the  continuity  equation 

(3-14)  W=d_mAVA)lA)r]' 

The  unclosed  conditional  expectation  (/(A,  VA)  |  A)  occurs  in  Eq.  (3.14)  because  the  derivative  VA 
couples  information  from  any  point  x  to  information  from  infinitesimally  nearby  points.  The  conditional 
expectation  projects  this  information  onto  local  information  at  x.  As  an  alternative,  it  is  possible  to  formulate 
the  continuity  equation  for  V  without  conditional  expectations,  but  only  as  a  coupled  equation  for  both  V 
and  the  two-point  joint  probability  density  V2(A,  A'),  where  A  =  A(x)  and  A'  —  A(x').  This  again  leads  to 
an  unclosed  hierarchy  of  equations  for  joint  PDF's  of  higher  order. 
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Conditional  expectations  do  not  appear  in  Craik’s  triad  model,  because  in  the  notation  of  the  model  Eq. 
(3.13),  /  depends  only  on  the  local  amplitude  A,  hence  trivially 

(3.15)  (f(A)  |  A)  =  /(A), 


which  allows  an  immediate  closure  for  the  Craik  triad  system.  On  the  other  hand,  assuming  that  the  matrix 
A  is  invertible,  the  evolution  equation  for  the  probability  density  P{A?)  corresponding  to  the  PSE  equations 
Eq.  (3.8)  is  given  by 


(3.16) 
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which  obviously  contains  unclosed  conditional  expectations. 

Here  we  will  only  outline  the  simplest,  most  basic  closure  scheme  for  the  PDF  equation  Eq.  (3.16),  the 
conditionally  Gaussian  closure  [13].  The  basis  of  the  closui'e  is  the  observation  that  in  a  Gaussian  random 
field,  conditional  expectations  containing  multipoint  quantities  are  all  closed  by  linear  regression,  which  is 
exact  for  Gaussian  fields.  To  illustrate  this  idea,  consider  the  simplest  case  of  a  zero  mean  homogeneous 
Gaussian  random  field  (f>(x).  The  Gaussian  property  implies  that  the  conditional  expectation  of  the  field 
value  (f>(x  +  r)  given  the  value  <j>{x)  satisfies  the  linear  regression  equation 


(3.17) 

Consequently, 

(3.18) 
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and  other  derivatives  can  be  evaluated  similarly.  The  conditionally  Gaussian  closure  simply  applies  results 
like  Eq.  (3.18)  to  a  general  random  field,  even  if  it  is  non-Gaussian.  This  step  can  be  compared  to  the 
quasi-normal  hypothesis  often  invoked  in  moment  closures. 

Of  course,  application  to  PSE  would  involve  inhomogeneous  random  fields.  The  requisite  modifications 
of  Eq.  (3.17)  are  presented  in  Appendix  II.  Here,  we  simply  quote  the  resulting  closure  for  the  PSE  equation, 


|  A_1D¥  +  A-1B  +  ^C(¥,  *)  '  W  *)[*  -  <*>] 


(3.19) 


&P  d_ 
dx  +  0 A? 

+A”1V  +  '  C(^’  ™ 
+A-1  (f(¥,  |^C(*,  ¥)  •  C(*,  ¥)[¥  -  <*)], 
|^C(*,  *)  C(*,  *)[*-<*>]* 


^1=0 


This  equation  contains  first  and  second  order  statistics  of  the  held  T'  which  are  defined  in  terms  of  the 
probability  density  V  by 


(3.20) 


(W)  =  j  d*  ’TP(’J') 
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and 


(3.21)  C(^,^)  =  J  d*  'P'FiP(’i') 

In  Eq.  (3.21),  the  product  is  an  outer  product. 

We  have  not  yet  addressed  the  issue  of  how  to  implement  an  integral  constraint  such  as  Eq.  (3.7)  in 
a  stochastic  formulation,  however,  further  comments  are  offered  in  Sect.  3.3  below.  Although  Eq.  (3.19) 
provides  a  complete  theory  of  the  evolution  of  the  joint  probability  density  function  of  modal  amplitudes,  it 
must  be  emphasized  that  it  is  based  on  the  assumption  of  the  conditionally  Gaussian  closure.  As  in  turbulence 
theory,  this  assumption  cannot  be  justified  a  priori ;  its  evaluation  can  only  be  based  on  comparisons  with 
experimental  or  numerical  data.  Although  not  be  discussed  in  this  work,  other  closure  schemes  based  on 
mapping  closure  or  on  Wierner-Hermite  expansions  (Y.  Kaneda,  private  communication)  are  also  possible. 
However,  they  will  all  require  a  significant  effort  in  terms  of  formulation,  validation,  and  refinement.  Thus, 
similar  to  the  triad  model  in  Sect.  2,  it  is  worth  investigating  the  application  of  the  particle  method  to  the 
PSE  system. 

3.3.  Stochastic  formulation:  truncation  of  modal  amplitude  equations  and  Langevin  models 
for  transition.  The  particle  method  will  require  calculating  a  large  number  of  individual  trajectories  with 
a  deterministic  model,  whether  the  PSE  or  other  deterministic  model.  Any  PSE  calculation  obviously  will 
require  a  truncation  of  the  Fourier  series  Eq.  (3.5)  to  a  finite  sum  of  modes.  However,  to  facilitate  the 
computation  of  the  required  number  of  particle  trajectories,  it  is  necessary  to  further  reduce  the  number 
of  modes  retained  in  whatever  deterministic  framework  is  applied  for  the  stochastic  analysis.  Accordingly, 
in  the  present  section,  we  consider  how  a  reduced  order  model  could  be  obtained  as  part  of  the  stochastic 
formulation.  This  amounts  to  deriving  ‘effective’  equations  governing  a  limited  number  of  explicitly  resolved 
modes.  Such  model  reduction  is  particularly  important  for  natural  transition,  which  would  contain  a  wide 
spectrum  of  modes  at  the  outset. 

For  simplicity,  consider  a  general  resonant  interaction  model 

(3.22)  Ai  —  jijAj  +  CijkAjAk  (i,j,  k  —  1,N) 

where  the  matrix  jij  is  diagonal,  so  that  there  is  no  linear  coupling  between  the  modes. 

Suppose  that  the  modes  can  be  divided  into  two  sets:  slow,  large  amplitude  resolved  modes  A+ ,  the 
modes  A*  with  1  <  i  <  Nr,  and  fast,  small  amplitude  unresolved  modes  A~,  the  modes  Ai  with  Nr  <  i  <  N . 
Let  it  be  required  to  model  the  evolution  of  the  complete  system  by  solving  modified  equations  for  the 
resolved  modes  alone.  It  is  possible  to  rewrite  Eq.  (3.22)  for  the  resolved  modes  in  the  form 

A^~  —  ('jip  +  CiprAr  )  A .1)7  A~  +  CirsAr  As 
(3-23)  —  JipAp  +  CipqApAq  +  fi 

where  we  have  introduced  the  effectively  modified  linear  growth  rates 

(3.24)  jfp  —  'yip  +  gip 
where 

(3.25)  gip  =  CiprA~ 
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and  the  forcing  attributed  to  unresolved  interactions 

(3.26)  fi  —  CirsA~A~ . 

At  this  point,  Eq.  (3.23)  is  nothing  more  than  a  formal  restatement  of  the  original  equations  Eq.  (3.22). 
But  since  the  unresolved  modes  are  assumed  to  vary  rapidly  compared  to  the  resolved  modes,  it  may  be 
reasonable  to  replace  the  terms  /;  and  gip  in  Eq.  (3.26)  by  random  quantities  with  statistics  to  be  evaluated 
from  the  resolved  part  of  the  motion.  Thus,  we  replace  the  term  by  a  random  force  with  correlation 
function 

(3.27)  ( fifj )  —  C irsC jim(Ar  As  Al  Am). 
and  the  quantity  yip  by  a  random  tensor  with  correlation  function 

(3.28)  ( 9ip9jq )  =  CiprCjqs(Ar  As  ) 

This  approximation  replaces  the  problem  of  explicitly  computing  the  evolution  of  the  unresolved  modes  with 
the  problem  of  modeling  of  the  correlations  on  the  right  sides  of  Eqs.  (3.27)  and  (3.28). 

Eq.  (3.23)  shows  that  the  unresolved  modes  modify  the  linear  growth  rates  and  act  as  a  random  force 
on  the  resolved  modes.  These  effects  will  reduce  the  phase  coherence  of  the  resolved  modes  and  could 
consequently  prevent  the  indefinite  growth  of  the  disturbance  amplitudes.  One  may  also  observe  this  on  the 
basis  of  the  PDF  evolution  equation  for  this  model.  In  the  absence  of  any  random  forcing,  the  PDF  evolution 
equation  only  contains  spatial  derivatives  of  the  convective  type  (see  Eq.  (2.12)).  Random  forcing  will  add 
diffusive  terms  to  this  equation  and  diffusion  will  counteract  the  tendency  toward  phase  coherence  exhibited, 
for  example,  in  Fig.  (2.4).  The  link  between  the  reduction  of  phase  coherence  and  the  suppresion  of  modal 
amplitude  growth  is  demonstrated  experimentally  in  Refs.  [9]  and  [29].  Another  viewpoint  on  this  link  in 
the  context  of  resonant  triad  theory  appears  in  Ref.  [48],  which  observes  that  higher-order  nonlinearity  of 
the  type  exhibited  in  Eq.  (3.2)  causes  amplitude-dependent  reduction  of  phase  coherence,  which  entirely 
suppresses  explosive  growth. 

To  investigate  the  connection  between  reduced  phase  coherence  and  inhibition  of  subharmonic  growth 
more  closely,  we  repeated  the  simpler  triad  computation  from  Ref.  [53]  after  modifying  the  phase  in  Eq.  (2.23) 
by  adding  white  noise  of  various  amplitudes  to  its  right  hand  side.  This  phase  randomization  is  equivalent 
to  randomization  of  the  coefficient  matrices  in  the  modified  PSE  equations  Eq.  (3.45).  Specifically,  we  set 

(3.29)  9  —  — — —  sin (cj)  —  9)  —  — — -  sin#  —  sin#  -p  Aw(x) 

Ti  r-2  r  i 

where  w(x )  is  a  white  noise  process  and  the  amplitude  A  was  set  equal  to  0.010,  0.050,  0.100,  0.200.  The 
degree  of  phase  decorrelation  caused  by  the  random  forcing  is  revealed  by  comparing  the  phase  evolution 
shown  in  the  left-hand  plots  in  Figs.  3. 1-3.4  with  the  phase  evolution  in  the  absence  of  random  forcing 
shown  in  Fig.  2.6.  At  the  smallest  noise  amplitude  A  —  0.010,  the  phase  locking  near  9  —  0  is  only 
weakly  perturbed:  the  phase  fluctuates  by  about  0.17T  about  the  deterministic  value.  At  this  level  of  phase 
randomization,  the  nonlinear  evolution  of  the  phase  dominates  the  effect  of  random  forcing.  But  when 
A  —  0.200,  it  is  apparent  that  the  phase  evolution  is  almost  entirely  dominated  by  the  random  forcing, 
and  the  phase  appears  to  evolve  as  a  Brownian  motion.  Correspondingly,  subharmonic  growth  is  almost 
entirely  suppressed  in  this  case.  Intermediate  values  of  the  noise  amplitude  delay  the  crossover  between  the 
amplitudes  of  the  subharmonics  and  the  primary. 
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We  stress  that  in  Figs.  3. 1-3.4,  the  suppression  of  energy  growth  in  the  subharmonics  has  been  ac¬ 
complished  entirely  through  phase  randomization:  energy  has  not  been  removed  from  the  system  directly 
in  these  calculations.  In  effect,  the  phase-locked  attractor,  which  promotes  subharmonic  growth,  has  been 
destroyed  by  stochasticity. 

Another  viewpoint  on  the  connection  between  phase  coherence  and  nonlinear  growth  is  through  the 
bicoherence  [29],  [35].  This  statistic  is  a  third-order  joint  moment  which  measures  the  correlation  between  one 
modal  amplitude  and  the  product  of  two  others.  Large  bicoherence  indicates  strong  triadic  coupling  between 
modes,  which  cannot  be  revealed  by  second-order  statistics.  Experimentally,  it  has  been  observed  that  in 
forced  transition,  the  bicoherence  spectrum  is  initially  strongly  peaked,  indicating  three-mode  resonance. 
As  transition  proceeds,  the  bicoherence  spectrum  broadens  and  weakens,  indicating  that  more  modes  are 
coupled  to  each  other,  but  that  the  degree  of  correlation  between  them  has  been  reduced.  In  future  work, 
we  will  consider  the  possibility  of  modeling  this  aspect  of  transition  through  stochastic  PSE  models,  which 
have  a  more  general  validity  than  the  present  stochastic  triad  model. 


Fig.  3.1.  Effect  of  random  noise  perturbation  of  phase  evolution  equation  on  subharmonic  amplification  during  resonant 
triad  evolution.  The  conditions  are  those  of  Fig.  2.6,  with  white-noise  perturbation  of  the  phase  evolution  equation  Eq.  (3.29 
with  amplitude  A  —  0.010. 

To  complete  the  reduced-order  modeling,  the  correlations  Eqs.  (3.27)  and  (3.28)  must  be  expressed  in 
terms  of  the  resolved  fields  alone.  As  an  example  of  a  possible  approximation,  suppose  that  the  amplitudes 
of  the  unresolved  modes  are  small  enough  that  nonlinear  interactions  among  them  can  be  ignored.  This 
approximation  leads  to  the  simple  “driven-mode”  approximation 

(3.30)  A-  =  t tfAj  +  CmAjA+. 

This  equation  can  be  solved  explicitly  in  the  form 

(3.31)  A~ ( x )  =  [  dy  Gar  (x  -  y)CVjk{y)Aj (y)A+  ( y ) 

Jo 

where 

(3.32)  Ga'{x-y)6we-***-*> 

is  the  linear  response.  The  substitution  of  Eq.  (3.31)  for  the  amplitudes  Aj  in  Eqs.  (3.27)  and  (3.28) 
achieves  the  goal  of  expressing  these  correlations  in  terms  of  the  resolved  modes  Af . 
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Fig.  3.2.  Effect  of  random  noise  perturbation  of  phase  evolution  equation  on  subharmonic  amplification  during  resonant 
triad  evolution.  The  conditions  are  those  of  Fig.  2.6,  with  white-noise  perturbation  of  the  phase  evolution  equation  Eq.  (3.29 
with  amplitude  A  —  0.050. 


Fig.  3.3.  Effect  of  random  noise  perturbation  of  phase  evolution  equation  on  subharmonic  amplification  during  resonant 
triad  evolution.  The  conditions  are  those  of  Fig.  2.6,  with  white-noise  perturbation  of  the  phase  evolution  equation  Eq.  (3.29 
with  amplitude  A  —  0.100. 


Because  it  uses  the  linear  response  function  Eq.  (3.32)  to  eliminate  the  unresolved  modes,  this  approx¬ 
imation  can  be  compared  to  the  quasi-normal  theory  of  turbulence.  A  refinement  is  offered  by  the  direct 
interaction  approximation  (DIA)  [23].  This  approximation  was  developed  as  a  statistical  closure  theory  for 
systems  with  a  large  number  of  nonlinearly  coupled  modes.  Applied  to  a  general  system  with  a  quadratic 
nonlinearity,  the  basic  assumption  behind  this  approximation  is  that  the  nonlinear  interactions  among  any 
specific  modal  triad  may  be  treated  as  small,  although  the  overall  effect  of  all  possible  triad  interactions  may 
nevertheless  be  quite  significant. 

The  DIA  improves  on  the  quasi-normal  theory  by  replacing  the  linear  response  function  of  Eq.  (3.32)  by 
a  generalized  response  function  defined  as  follows.  First,  evaluate  the  linearized  equation  for  the  response  of 
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Fig.  3.4.  Effect  of  random  noise  perturbation  of  phase  evolution  equation  on  subharmonic  amplification  during  resonant 
triad  evolution.  The  conditions  are  those  of  Fig.  2.6,  with  white-noise  perturbation  of  the  phase  evolution  equation  Eq.  (3.29 
with  amplitude  A  —  0.200. 


A  to  an  infinitesimal  perturbation  £*  as 

(3.33)  5A^  —7 ijSAj  A  CijkdAj  Af,  +  CijkA-  5A^  +£$. 

The  generalized  response  function  is  defined  as  the  average  ‘sensitivity  matrix’ 

lx  4- 

(3.3d)  G'7  -  ^ 

which  replaces  the  linear  response  function  Eq.  (3.32)  in  Eq.  (3.28).  The  theory  is  completed  by  an  evolution 
equation  for  the  response  matrix  Gij.  In  the  context  of  transition  theory,  the  generalized  response  function 
can  perhaps  be  compared  to  the  quantities  a  of  the  PSE:  unlike  linear  theories  based  on  the  Orr-Sommerfeld 
equation,  PSE  determines  these  functions  in  a  fully  nonlinear,  coupled  fashion.  Similarly,  the  DIA  response 
function  replaces  the  linear  response  by  a  quantity  which  accounts  for  nonlinear  interactions. 

Note  that  the  integration  in  Eq.  (3.32)  introduces  history  effects,  similar  to  the  high-Reynolds  number 
asymptotic  models.  This  again  introduces  closure  difficulties  for  stochastic  formulations  of  the  amplitude 
evolution  equations. 

3.4.  Link  to  turbulence  modeling.  The  results  of  Sects.  3.2  and  3.3  will  now  be  applied  to  the  PSE 
equations.  Our  goal  is  to  develop  a  usable  stochastic  form  of  the  PSE  equations  in  which  only  a  reasonably 
small  set  of  modes  is  resolved  explicitly,  and  the  effects  of  the  remaining  modes  is  modeled  through  random 
forces  and  random  coefficients.  This  model  will  provide  a  natural  interface  between  transition  modeling  and 
turbulence  modeling  if  the  unresolved  motion  can  be  characterized  statistically  by  a  turbulence  model. 

The  breakdown  of  PSE  calculations  shortly  after  the  rise  in  the  skin  friction  coefficient  has  been  discussed 
in  Sect.  3.2.1.  To  continue  the  PSE  past  this  point,  and  perhaps  even  into  the  fully  turbulent  region,  we 
will  assume  that  the  PSE  model  would  be  adequate  if  the  basis  of  modes  could  be  expanded  indefinitely; 
but  rather  than  explicitly  resolving  the  new  modes  generated  during  the  phase  of  rapid  spectral  broadening 
by  PSE,  we  will  treat  them  as  unresolved  components  of  the  motion  to  be  modeled  following  Eq.  (3.23). 
Since  these  modes  are  both  faster  in  time  scale  and  smaller  in  amplitude  than  the  resolved  modes,  this  type 
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of  modeling  is  reasonable.  Our  contention  will  be  that  a  randomized  PSE  model  will  provide  two  required 
parts  of  an  integrated  turbulence  and  transition  model,  namely  a  transition-like  model  which  applies  into 
the  breakdown  region,  and  which  interfaces  naturally  with  a  RANS  model  which  will  be  applied  in  the  fully 
turbulent  region. 

Following  the  notation  of  Sect.  3.3,  we  decompose  the  disturbance  into  resolved  and  unresolved  modes 
as 

fi=  E  A  tn^,y)ei{m0z~rwt)  +  E 

\m\<M,\n\<N  \m\>M,\n\<N 

(3.35)  +  E  A-Jx,y)eiW>-™t)+  E  A^Jx^e^^ 

\m\<M\\n\>N  \m\>M  ,\n\>N 

where  the  +  modes  will  be  resolved  explicitly  and  the  —  modes  will  be  modeled.  After  applying  the  wave- 
amplitude  decomposition  Eq.  (3.5),  the  resolved  modes  are  found  to  satisfy  the  equations 


d*+ 


hIL  -  w 

—  v  n 


d2V+ 


^  +  fro,„(*+  V*+)  +  Fra,n(¥-,  V¥") 


(3.36) DTOjn^,^  rt  +  Amjn  ^  I  ~m,n  o  —  -  m,n  o  2 

where  the  contributions  from  the  unresolved  modes  appear  in  the  term  FTOj„  which  is  defined  by 

(3.37)  Fm,n(^-,v®-)  -  v$+)  vr)  +  W*-,  Vf) 

and  the  terms  on  the  right  side  of  Eq.  (3.37)  are  defined  by 

(3.38) 


ij  ,m,n 


(3.39) 

(3.40) 

where 

(3.41) 


J  ,m,n 


f^n(*-,V¥+)  =  E 

i+ir  —m,j+jr  —n 

f ™,n(®+>V®")=  E 

—m,j+jr  —n 

fra,n(1®r  ?  =  ^  ^  ir ,jr 

i+ir  —n 

,jr  iTTiiTi  —  GXp  %  I  )  +  Oi'fix  )  Ot m,n(, &  'jdx 

.  Jxq 


For  brevity  of  notation,  restrictions  on  the  summation  ranges  in  Eq.  (3.40)  which  arise  from  the  definitions 
in  Eq.  (3.35)  have  been  omitted  in  Eqs.  (3.38)-(3.40). 

The  modified  PSE  Eq.  (3.37)  can  be  rewritten  as 


d^f+r  , 

rrv  ,nr 


(3.42) 
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E  D'  vEr+  -j-  A'  m',n'  4-  B' 
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dy~ 


where  the  new  convective  coefficients  include  the  effects  of  unresolved  modes: 


(3.43) 

(3.44) 

(3.45) 
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and  the  forcing  term 


(3.46) 


Qm,n  —  f/n , u (tit  ,  VW  ) 


contains  contributions  from  unresolved  modes  alone. 

Eqs.  (3.42)-(3.46)  simply  restate  the  PSE  equations  with  a  larger  set  of  modes.  The  goal  now  is  to 
model  the  effect  of  the  added  fields  dt-  without  explicitly  resolving  them,  as  in  large  eddy  simulations  of 
turbulence.  One  approach  would  be  to  develop  approximate  equations  of  evolution  for  the  unresolved  scales 
in  terms  of  the  resolved  scales.  Instead,  it  is  also  possible  to  treat  the  random  couplings  in  Eqs.  (3.43)~(3.45) 
and  the  term  qTOi„  in  Eq.  (3.46)  as  Gaussian  random  processes  with  statistics  determined  by  the  resolved 
fields  as  discussed  in  Sect.  3.2.  When  this  is  done,  the  original  PSE  is  replaced  by  a  Langevin-like  model 
which  includes  a  random  force.  Unlike  a  strict  Langevin  model,  in  which  a  random  force  is  added  to  a 
deterministic  equation,  the  model  proposed  here  will  have  both  a  random  force  and  random  coefficients. 
Applications  of  Langevin  models  to  turbulence  theory  have  been  proposed  by  [31]. 

The  random  force  qm>„  has  the  correlation  function 

(3.47)  (qro.n  q— rra,— «.}  —  Vd^j,]^ 

h+h' =— m,  l-\-V=—n 

Note,  first,  that  Eq.  (3.47)  defines  a  space-time  correlation  and,  second,  the  process  q  is  not  white  noise  in 
time. 

We  will  outline  one  particular  model,  suggested  by  Sect.  3.2,  to  close  the  PSE  equations  for  resolved 
modes.  To  begin,  replace  the  complete  PSE  equation  for  the  unresolved  modes,  following  Eq.  (3.23)  by 


(3.48) 


^m,n^ m,n  T  An 


<9d>~ 

m>n  i  p 
i  o  I  &n 
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i-\-ir  —m,j+jr  —n 


In  making  this  approximation,  we  argue  that  the  disturbances  represented  by  the  mode  amplitudes  \P~  are 
small  compared  to  the  resolved  disturbances  dt.  Then  it  is  reasonable  to  ignore  the  nonlinearity  of  the  d* 
evolution  equations  and  treat  the  dW  modes  as  driven  by  the  resolved  modes. 

The  matrices  A  •  •  •  D  in  Eq.  (3.48)  contain  the  complex  growth  rates  aTO,n  which  are  found  in  PSE  as 
part  of  the  modal  evolution.  Since  the  dt  modes  are  not  to  be  explicitly  resolved,  one  might  replace  these 
quantities  by  the  linear  growth  rates  obtained  by  linearizing  the  governing  equations  about  the  resolved 
disturbance 


(3.49)  u°  =  A 

\m\<M\\n\<N 

With  this  prescription  for  the  complex  growth  rates,  Eq.  (3.48)  is  determinate,  and  the  force  correlation  Eq. 
(3.47)  can  be  evaluated  in  terms  of  the  resolved  disturbance  alone. 

By  specializing  to  the  (m,  n)  —  (0, 0)  mode,  we  can  form  the  equation  for  the  mean  flow  modification.  In 
this  equation,  the  amplitude  of  the  random  perturbation  d/^y  can  be  expected  to  grow  in  comparison  to  the 
resolved  stresses  qo,o-  As  noted  earlier,  the  sum  \Py  0  +  q0jo  defines  the  Reynolds  stress  gradients.  Although 
the  PSE  Langevin  model  is  a  completely  consistent  turbulence  theory,  in  the  interests  of  economical  calcu¬ 
lations,  it  can  be  replaced  by  a  turbulence  model  once  the  random  contribution  dominates  the  contribution 
from  the  resolved  modes.  Relevant  switchover  strategies  of  this  type  have  been  proposed  in  a  somewhat 
different  context  by  Spalart  [42]  for  detached  eddy  simulation  and  by  Speziale  [43]  for  combining  DNS  with 
RANS. 
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It  is  also  of  some  interest  to  note  that  the  PSE  equations  include  a  kind  of  two-equation  model.  To  derive 
this  model,  we  replace  the  modal  amplitude  vector  \P  by  a  vector  tP  in  which  the  amplitude  components 
pertain  only  to  the  velocities.  Indeed,  multiplying  Eq.  (3.8)  by  gives  a  disturbance  energy  equation, 
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Since  \P  is  a  vector  and  the  quantities  A  •  •  •  D  are  matrices,  the  products  in  Eq.  (3.50)  are  all  scalars.  The 
effective  dissipation  rate  in  this  equation  is  the  last  term,  its  equation,  also  found  from  Eq.  (3.8) 
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Note  that  the  above  dissipation  rate  represents  only  the  nonlinear  energy  transfer  from  large  to  small  scales. 
It  ignores  viscous  dissipation;  in  that  sense,  it  may  differ  from  the  total  dissipation  rate  used  in  conventional 
two-equation  models. 

We  conclude  by  briefly  discussing  how  the  PDF  equation  for  PSE  modal  amplitudes  can  be  modified  to 
account  for  the  random  force  introduced  above  as  a  model  for  the  unresolved  modes.  A  complete  PDF  model 
for  the  stochastic  PSE  cannot  be  formulated  exactly,  because  of  the  implicit  procedure  used  to  determine  the 
growth  rates  a.  The  evolution  of  these  quantities  is  highly  path  dependent,  and  it  is  impossible  to  construct 
a  stochastic  model  of  the  evolution  of  all  possible  transition  paths.  Also,  the  appearance  of  stochastic 
coefficients  will  require  extensions  to  the  standard  methods.  Thus,  the  present  discussion  will  be  limited  to 
a  very  special  case,  namely  we  will  assume  that  the  a  can  be  determined  once  and  for  all  for  all  paths,  and 
only  the  effect  of  the  random  force  term  q  will  be  considered. 

Using  the  conditionally  Gaussian  closure,  the  modified  equations  will  have  the  form 
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Again,  the  crucial  effect  of  the  random  force  appears  in  the  last  two  lines  of  Eq.  (3.52)  which  are 
diffusive  terms  indicative  of  the  weakening  of  modal  correlations.  These  diffusive  terms  again  have  the 
effect  of  preventing  the  indefinite  buildup  of  modal  phase  correlations,  a  critical  feature  of  any  theory  which 
proposes  to  integrate  the  dynamics  of  transition  and  of  turbulence.  We  stress,  however,  that  numerical 
solution  of  Ecp  (3.52)  is  probably  impractical. 

Could  a  model  like  Eq.  (3.52)  be  used  in  practice?  The  large  number  of  degrees  of  freedom  could  perhaps 
be  handled  by  the  particle  method,  by  finding  a  simpler  system  than  the  stochastic  PSE  equations  for  which 
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the  PDF  of  the  modal  amplitudes  is  also  governed  by  Eq.  (3.52).  This  recommendation  of  Pope  [31]  has 
been  succesfully  applied  to  problems  like  reacting  hows,  which  in  a  Lagrangian  frame  reduce  to  ordinary 
differential  equations.  In  the  present  problem,  the  ‘paths’  are  evolving  fields  described  by  partial  differential 
equations.  Extension  of  the  particle  method  to  this  class  of  problems  will  require  further  generalization  of 
the  basic  particle  method  formulation.  The  more  fundamental  problem  for  stochastic  formulation  of  the 
PSE  is  that  each  transition  path  will  correspond  to  a  distinct  evolution  of  the  complex  wavenumbers  am>n. 
Thus,  even  more  approximations  will  be  needed  to  formulate  the  PSE  in  a  fully  stochastic  setting. 

4.  Summary  and  future  work.  Integrated  modeling  of  transition  and  turbulence  has  emerged  as  a 
crucial  need  in  many  practical  aerodynamic  applications.  To  retain  the  physical  basis  necessary  to  explain 
the  relevant  parametric  dependencies  and,  hence,  achieve  the  desired  prediction  accuracy,  such  integrated 
computations  must  be  done  seamlessly  over  the  entire  transition  process,  from  the  initial  disturbances  to  fully 
developed  turbulence.  Thus,  two  levels  of  integration  are  inherent  to  an  approach  of  this  type,  a  coupled 
prediction  of  the  laminar-through-transitional  and  turbulent  regions  of  the  flowfield  (at  the  top)  and  an 
integrated  prediction  of  the  various  stages  of  transition  (on  the  tier  underneath).  In  addition,  it  is  desirable 
to  pose  the  transition  prediction  problem  in  a  stochastic  context,  in  order  to  both  facilitate  an  easy  interface 
with  the  statistical  models  for  turbulence  and  to  quantify  the  effect  of  randomness  and/or  uncertainties 
associated  with  the  input  quantities  that  are  required  for  transition  prediction.  The  present  work  may  be 
regarded  as  an  embryonic  attempt  to  explore  the  above  needs,  with  an  emphasis  on  the  stochastic  aspects  of 
nonlinear  disturbance  evolution  in  a  laminar  boundary  layer  and  its  potential  continuation  into  the  regions 
of  laminar  breakdown  as  well  as  the  fully  turbulent  flow  downstream. 

The  fundamental  properties  of  resonant  wave  interactions  in  boundary  layer  flows  are  that  subharmonics 
can  grow  rapidly  through  parametric  excitation,  and  that  parametric  growth  is  typically  followed  by  an 
explosive  nonlinear  growth  of  all  modal  amplitudes.  The  explosive  growth  phase  can  be  delayed  if  the  initial 
phases  are  suitably  related.  Our  calculations  confirm  these  findings  by  previous  researchers  and  quantify 
the  dependence  of  disturbance  amplitudes  on  the  initial  conditions  by  evaluating  the  evolution  of  their  joint 
probability  density  function  for  a  given  description  of  the  initial  disturbances  not  in  deterministic  terms,  but 
in  terms  of  corresponding  statistical  measures. 

A  basic  objection  to  both  linear  stability  theory  and  resonant  interaction  theory  is  that  the  predicted 
indefinite  amplification  of  disturbances  is  not  observed,  even  though  the  mean  flow  could  in  principle  provide 
an  energy  source  for  such  growth.  ‘Nonlinearity,’  without  further  qualification,  is  sometimes  cited  as  the 
reason  for  the  saturation  of  the  disturbance  amplitude;  but  even  the  spreading  of  the  disturbance  over  every 
more  modes  is  not  an  explanation.  Instead,  the  prevention  of  the  development  of  strong  phase  correlations 
by  coupling  to  a  larger  set  of  modes  appears  to  have  a  significant  role  as  well.  Therefore,  in  developing  an 
integrated  theory  of  transition  and  turbulence,  it  will  be  important  to  recognize  the  proper  role  of  phase 
correlations  in  determining  the  onset  of  transition. 

The  probabilistic  analysis  has  been  extended  to  the  PSE  equations,  which  are  to  date  the  most  compre¬ 
hensive,  albeit  non-rigorous,  model  for  engineering  prediction  of  laminar-turbulent  transition.  The  crossover 
from  the  PSE  equations  to  a  turbulence  model  has  been  described  through  a  PSE  Langevin  model,  in  which 
the  PSE  equations  are  modified  by  random  forces  and  random  dephasing  terms.  When  the  random  per¬ 
turbations  grow  sufficiently,  the  PSE  Langevin  model  can  be  replaced  by  a  RANS  model  in  the  interest  of 
economy  of  computation. 

Our  conclusion  is  that  a  stochastic  formulation  of  the  PSE  transition  model  can  link  transition  modeling 
and  turbulence  modeling  in  the  required  seamless  fashion.  Stochastic  closures  for  unresolved  motions  are  the 
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most  important  requirement.  The  recommendations  in  Sect.  3.3  were  based  directly  on  transition  theory, 
but  an  interesting  alternative  (which  may  be  easily  implemented  in  practice)  would  be  to  apply  the  ideas  of 
large  eddy  simulation  to  formulate  analogs  of  the  Smagorinsky  and  dynamic  Smagorinsky  model  appropriate 
to  transitional  flows  within  the  PSE  framework.  This  approach  has  the  potential  to  make  the  crossover  to 
turbulence  models  much  more  straightforward. 

The  primary  motivation  behind  considering  reduced-order  models  for  PSE  was  to  facilitate  engineering 
predictions  through  the  transition  region.  However,  another  benefit  would  be  to  provide  detailed  physical 
insights  into  the  transition  process  via  a  coupled  PSE-LES/DNS  approach.  Essentially  by  extending  the 
PSE  predictions  further  into  the  laminar  breakdown  region,  the  computational  domain  for  LES/DNS  could 
significantly  shortened  (see,  for  example,  [19],  [32]).  Indeed,  if  PSE  calculations  can  be  successfully  continued 
even  into  the  turbulent  region,  there  is  the  interesting  possibility  of  being  able  to  predict  sublayer  dynamics 
including  wall  streaks  and  bursting  frequencies.  Application  of  PSE  in  this  respect  offers  an  attractive 
alternative  to  other  modal  analyses  of  fully  turbulent  flows  which  have  been  pursued  over  many  years  ([54], 
[47]).  Of  course,  to  make  this  application  feasible,  extensive  comparisons  against  numerical  and,  where 
available,  experimental  databases  will  have  to  be  earned  out.  We  also  plan  to  conduct  an  a  priori  analysis  of 
the  simulation  data,  both  in  order  to  evaluate  the  accuracy  of  candidate  models  for  residual  stress  components 
and  to  ascertain  the  validity  of  PSE  methodology  past  the  onset  of  transition. 

So  far,  we  have  focused  our  attention  on  stochastic,  integrated  prediction  of  the  transition  process.  The 
next  phase  of  this  work  will  involve  linking  transition  prediction  and  conventional  CFD  codes.  In  brief,  the 
overall  computation  involves  an  iterative  procedure  consisting  of  a  RANS  computation  over  the  entire  flow 
field  coupled  with  an  ‘unsteady’  computation  (PSE  and  residual  stress)  over  the  laminar  and  transitional 
region.  The  mean  flow  (obtained  directely  from  RANS  or  via  a  boundary  layer  solver)  is  fed  from  the  RANS 
model  to  the  transitional  module,  which  predicts  the  extent  of  the  laminar  and  transitional  regions,  plus  the 
mean  Reynolds  stress  distribution.  This  information,  in  turn,  is  injected  into  the  RANS  simulation.  This 
procedure  is  analogous  to  ‘zonal’  modeling  for  fully  turbulent  flow  fields,  except  that  the  zonal  boundaries 
are  not  predetermined  in  our  case  (also,  the  unsteady  caleualtion  may  take  place  over  a  different  grid  using 
spatial  marching  approach  instead  of  a  pseudo-time  marching  typical  of  RANS  calculations).  A  robust  im¬ 
plementation  of  this  procedure  will  require  special  attention  to  numerical  convergence  of  the  global  iteration 
procedure.  Other  related  numerical  issues  have  recently  been  investigated  in  [44]. 

5.  Acknowledgments.  The  authors  gratefully  acknowledge  helpful  technical  discussions  with  Drs.  C. 
L.  Streett  and  C.  L.  Chang. 
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6.  Appendix  I:  Finite  volume  integration  of  PDF  evolution  equation.  In  this  Appendix,  we 
present  a  brief  description  of  the  numerical  method  used  to  solve  the  PDF  evolution  equation  Eq.  (2.12)  for 
a  resonant  triad. 

We  apply  finite  volume  discretization,  dividing  the  space  of  the  independent  variables  r,  r3 ,  9  into  rect¬ 
angular  cells  with  side  lengths  Ar,  A r3 ,  A 9,  where 

(6.1)  N(x  |  x0,a)  =  —l=e-(x-xo)2/^2 

C'2ira~ 


and  assuming  that  P  is  constant  in  each  cell.  Let  cell  ( i,j,k )  be  bounded  by  the  planes  r  —  (i  +  1)A r,r  — 
iAr,  r-s  —  ( j  +  l)Ar3,  r3  —  jAr-j,  9  —  (fc  +  1)A0,  9  —  kAO  Let  Pijtk  denote  the  value  of  P  in  this  cell,  and  let 
7 k  denote  the  value  of  P  on  the  face  on  which  r  —  iAr,  with  analogous  definitions  for  the  values  7r3  and 
.  Integrating  Eq.  (2.12)  over  each  cell  produces  the  discrete  equation 


(0.2) 


Pi,j,k  if  +  At)  —  Pi,j,kif)  "F 


At 


A 14 


i,j,k 


{Fi,j,k  + 


where 


-  /  dr  | :(FrP ) 

=  Fr((i  +  l)Ar,  j A7'3,  kA9)-Kri+1j  k  -  Fr(iAr,jAr3,kA0)irriJik 

Fh,k  =  J  dr3  -^{FzP) 

=  F3(iAr,  ( j  +  l)Ar3,  kA9)  F3{iAr,  j  Ar  3,  kA9)Trfjk 

Fij,k  =  Jd9-^(F0P) 

(6.3)  =  Fff(iAr,jAr3,  ( k  +  l)A9)'K(,i  j  k+1  -  Fe(iAr,jAr3,  kA9)-K9i  j  k. 


The  assumption  that  P  is  constant  in  each  cell  implies  that  the  fluxes  T  can  be  evaluated  exactly,  but 
the  problem  remains  to  evaluate  the  values  on  the  faces  it  by  interpolating  the  values  of  P  in  each  cell.  The 
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most  obvious  procedure  of  averaging  the  values  of  P  across  each  face: 


(6.4) 


wi,j,k  —  2  Pi—1  ,j,k\ 

3  1 

P,j,k  ~  +  Pi,j—l,k ] 

wi,j,k  —  2  hk  Piiiik— l] 


is  well-known  to  be  unstable  [17].  A  stable  scheme  results  from  the  “upwind”  method 


(6.5) 


7T. 


i,j,k 


P(i,j,k)  if  Fr(i,j,  k)  >  0 
P{i  -  l,j,  k)  if  Fr(i,j,  k )  <  0 


with  analogous  definitions  for  7 r3  and  ,  in  which  the  value  of  ir  is  the  value  of  P  at  the  cell  into  which  the 
corresponding  flux  P  points. 

The  accuracy  of  the  computation  degrades  quickly  if  there  is  a  probability  flux  across  the  boundaries 
of  the  region  of  integration.  Since  the  amplitude  distributions  spread  out,  the  region  of  integration  must 
increase  to  prevent  any  loss  of  total  probability.  To  maintain  the  total  probability,  the  box  sizes  are  rescaled 
at  each  time-step  so  that  a  fixed  number  of  boxes  always  contain  all  but  a  small  specified  amount  of  the 
total  probability.  The  rescaling  is  done  on  the  basis  of  the  standard  deviations  of  the  marginal  distributions 
P(r)  and  P(r3).  Strictly  speaking,  rescaling  based  on  the  standard  deviation  is  accurate  only  for  Gaussian 
distributions,  for  which  the  region  from  about  —4.5(7  to  +4.5(7  contains  almost  all  of  the  probability;  however, 
the  calculations  showed  that  this  simple  procedure  was  satisfactory  in  practice. 

To  implement  this  variable  box-size,  the  governing  equation  is  modified  as  follows:  define  the  rescaled 
variables 


r  —  crrx  +  TO 

(6.6)  r-i  -  (t-sx-s  +  to3. 

Note  that  since  the  problem  is  periodic  in  9,  rescaling  of  9  is  not  possible.  Instead  of  Eq.  (2.12)  we  solve  the 
modified  equation 

<6'7)  W  +  ([w  ,h)p])  +  +  ^  (l(F»  - &3X3  - ,h3}Pl)  +  re(F,p)  =  0 

where  the  fluxes  Fr ,Fs,Fo  are  the  same  as  the  fluxes  in  Eq.  (2.6)  and 


m  — 

J 

f  FrPdV 

7h 3  =  j 

[  F'iPdV 

&  — 

J 

f  xFrPdV 

(6.8) 

&3  =  j 

f  x3F3PdV 

For  parametric  excitation,  the  governing  equation  Eq.  (2.12)  is  modified  by  altering  the  flux  terms  to 

Fr  —  rr.i  cos  9  +  or 
Fi  —  <rr3 

(6.9)  Fff  —  —  2r3sin# 

with  the  second  modification  Eq.  (6.8)  for  rescaling,  but  otherwise,  the  numerical  integration  method  is 
unchanged. 
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7.  Appendix  II:  Results  for  inhomogeneous  Gaussian  random  fields.  In  this  Appendix,  we 
derive  some  results  used  to  formulate  the  conditionally  Gaussian  closure  for  PDF  evolution  equations. 

For  an  inhomogeneous  Gaussian  random  vector  field  4>i(x),  introduce  the  statistics 

ai(x)  = 

rij(y,x)  =  (tpi(y)ipj(x)} 

Q j(y,x)  =  {fi(y)){fj{x)}  -  {ipi(y)%l>j (x)} 

(7.1)  =  -(U>j(y)  -  W’j{yMi’i(x)  -  {i>i(x)}]}. 

Note  that  the  covariance  matrix  C  is  nonsingular. 

Conditional  expectations  are  found  from  the  lineal-  regression  equation  which  is  exact  for  Gaussian 
random  fields 


(7.2)  (il>i(y)  \  ■  ■  -il>n(x)}  =  a,i(y)  +  Qim(y,x)Qmj{x,x)  1[4’j{x)  -  a,j{x)] 


Consequently,  conditional  expectations  of  derivatives  are  given  by 

I  Vhfc) '  ■■4’n(x)^  =  a'i(x)  +  ^^-( y ,x)  | y=x  Qmj{x,x)-l[4>j{x)  -  dj{x)] 
(7-3)  |  ■  ••4’n(x)j  =  a"{x)  +  ^7  (V-.  x)  \y=x  C  mj(x,x)-1[i(>j(x)  -  <ij(x)] 


where  we  can  substitute 


(7.4) 


d 2 


dy~  ~~lJ 


Qj(y,x)  \y=x  — 


(a,'{x))(aj(x)} 


d_ 

dx 


d4’i{x)  \  f  &4>i{x)  dil>j{x)\ 

dx  Vj(  ])  \  dx  dx  / 


so  that  multiplication  by  the  viscous  term  V  leads  to  a  dissipative  correlation. 


34 


